$title sage is an applied general equilibrium model



* ------------------------------------------------------------------------------
* set options
* ------------------------------------------------------------------------------

* set default benchmark file
$if not set benchmark_file $set benchmark_file "data\default_aggregation.gdx"

* set default for disaggregating activities based on extant capital
$if not set putty_clay $set putty_clay 1

* set default for parameter file
$if not set parameter_file $set parameter_file "data\parameters.gms"

* by default don't use balanced growth starting values if a baseline is supplied
$if not set balanced_start_values $set balanced_start_values 0

* optional additive perturbation of the y(t,r,s) starting values
$if not set perturb_start $set perturb_start 0

* zero out the baseline file if one is not provided
$if not set gdx_baseline_file $set gdx_baseline_file 0

* zero out the policy file specification if one is not provided
$if not set policy_file $set policy_file 0

* save the results in a gdx file by default
$if not set gdx_save $set gdx_save 1

* default name for the gdx file where the results are saved
$if not set gdx_results_file $set gdx_results_file "output\results.gdx"

* set default name for output file if needed
$if not set output_file $set output_file "output\results.csv"

* zero out the extra code file specification if one is not provided
$if not set extra_code_file $set extra_code_file 0

* the following reduce the size of the listing file
$offsymxref
$offsymlist
option limcol = 0;
option limrow = 0;
option solprint = off;



* ------------------------------------------------------------------------------
* define model sets
* ------------------------------------------------------------------------------

* basic model sets
sets
  s                     sectors
  r                     regions
  h                     households
  trd                   trade markets
  f                     primary factors
  t                     time periods
  tfirst(t)             first time period
  tlast(t)              last time period;

* subsets
sets
  ftrd(trd)             foreign trade market
  dtrd(trd)             national trade market
  sm(s)                 materials sectors
  sen(s)                primary energy sectors
  sel(s)                electricity sector
  sstd(s)               standard production sectors
  sres(s)               resource extraction sectors
  sagf(s)               agriculture and forestry sectors
  scm(s)                consumption materials sectors
  stn(s)                transportation sector
  scen(s)               consumer primary energy sectors;

* set aliases
alias
  (s,ss), (r,rr), (h,hh), (t,tt);



* ------------------------------------------------------------------------------
* define model parameters
* ------------------------------------------------------------------------------

parameters

* benchmark data
  y_ex0(r,s)            benchmark - output with extant capital
  ld_ex0(r,s)           benchmark - labor demand for use with extant capital
  kd_ex0(r,s)           benchmark - extant capital demand
  id_ex0(r,ss,s)        benchmark - intermediate demand for use with extant capital
  res_ex0(r,s)          benchmark - resource for resource extraction sectors with extant capital
  y0(r,s)               benchmark - output with new capital
  klem0(t,r,s)          benchmark - capital(k)-labor(l)-energy(e)-materials(m) bundle
  kle0(t,r,s)           benchmark - capital(k)-labor(l)-energy(e)-materials(m) bundle
  kl0(t,r,s)            benchmark - capital(k)-labor(l) bundle
  ene0(t,r,s)           benchmark - primary energy(en)-electricity(e) bundle
  en0(t,r,s)            benchmark - primary energy(en) bundle
  mat0(t,r,s)           benchmark - materials(m) bundle
  ld0(r,s)              benchmark - labor demand
  kd0(r,s)              benchmark - capital demand
  id0(r,ss,s)           benchmark - intermediate demand
  ld_base(t,r,s)        benchmark - labor demand
  kd_base(t,r,s)        benchmark - capital demand
  id_base(t,r,ss,s)     benchmark - intermediate demand
  res0(r,s)             benchmark - fixed factor resource
  a0(r,s)               benchmark - armington aggregate
  dn0(r,s)              benchmark - domestic(d)-national(n) bundle
  d0(r,s)               benchmark - domestic supply
  m0(r,s,trd)           benchmark - imports
  x0(r,s,trd)           benchmark - exports
  cl0(r,h)              benchmark - consumption(c)-leisure(l) bundle
  leis0(r,h)            benchmark - leisure
  c0(r,h)               benchmark - consumption bundle
  cd0(r,s,h)            benchmark - consumer demand
  cd_base(t,r,s,h)      benchmark - consumer demand
  cm0(t,r,h)            benchmark - consumption materials bundle
  cene0(t,r,h)          benchmark - consumption electricity and energy bundle
  cen0(t,r,h)           benchmark - consumption primary energy bundle
  cme0(t,r,h)           benchmark - consumption materials and energy bundle
  inv0(r)               benchmark - investment bundle
  i0(r,s)               benchmark - investment demand
  invh0(r,h)            benchmark - household investment
  gov0(r)               benchmark - government bundle
  g0(r,s)               benchmark - government demand
  g_base(t,r,s)         benchmark - government demand
  le0(r,h)              benchmark - labor rental
  l0(r,h)               benchmark - household labor supply
  re_ex0(r,h)           benchmark - extant capital rental income
  re0(r,h)              benchmark - new capital rental income
  rese_ex0(r,s,h)       benchmark - resource income from extant capital
  rese0(r,s,h)          benchmark - resource income from new capital
  r0(r)                 benchmark - new capital rental
  bopdef(r,h)           benchmark - international balance of payments
  incadj0(r,h)          benchmark - government to household lump sum transfer
  res_share(r,s)        benchmark - fixed factor share of capital
  n0(r,h)               benchmark - households

* benchmark tax rates
  tl0(r)                tax rate - labor use
  tk0(r)                tax rate - capital use
  tc0(r)                tax rate - consumption
  ty_ex0(r,s)           tax rate - production using extant capital
  ty0(r,s)              tax rate - production using new capital

* tax rates
  tl(t,r)               tax rate - labor use
  tk(t,r)               tax rate - capital use
  tc(t,r)               tax rate - consumption
  ty_ex(t,r,s)          tax rate - production using extant capital
  ty(t,r,s)             tax rate - production using new capital

* elasticities
  se_klem(s)            substitution elasticity - materials and value-added-energy bundle
  se_kl(s)              substitution elasticity - value-added
  se_kle(s)             substitution elasticity - value-added-energy bundle
  se_ene(s)             substitution elasticity - electricity and primary energy bundle
  se_en(s)              substitution elasticity - primary energy bundle
  se_rklem(t,r,s)       substitution elasticity - r-klem bundle in resource extraction sectors
  se_dn(s)              substitution elasticity - domestic-national armington elasticity
  se_nf(s)              substitution elasticity - national-foreign armington elasticity
  se_cl(r,h)            substitution elasticity - consumption-leisure bundle
  se_c                  substitution elasticity - consumption bundle
  se_cme                substitution elasticity - consumption materials and energy bundle
  se_cen                substitution elasticity - consumption electricity and energy bundle
  se_cene               substitution elasticity - consumption primary energy bundle
  se_cm                 substitution elasticity - consumption materials bundle
  eta(r,h)              substitution elasticity - inverse of consumption over time
  te_dx(s)              transformation elasticity - domestic and exported production
  te_k_ex               transformation elasticity - sector differentiated extant capital

* benchmark cost shares
  cs_y_ex(r,*,s)        cost share - input shares for use with extant capital
  cs_rklem(t,r,*,s)     cost share - resource(r)-capital(k)-labor(l)-energy(e)-materials(m) bundle
  cs_klem(t,r,*,s)      cost share - capital(k)-labor(l)-energy(e)-materials(m) bundle
  cs_kle(t,r,*,s)       cost share - capital(k)-labor(l)-energy(e) bundle
  cs_kl(t,r,*,s)        cost share - capital(k)-labor(l) bundle
  cs_ene(t,r,*,s)       cost share - primary energy(en)-electricity(e) bundle
  cs_en(t,r,*,s)        cost share - primary energy(en) bundle
  cs_mat(t,r,sm,s)      cost share - materials(m) bundle
  cs_dn(r,*,s)          cost share - domestic(d)-national(n) armington aggregate
  cs_nf(r,*,s)          cost share - national(n)-foreign(f) armington aggregate
  cs_dx(r,*,s)          cost share - destination differentiated output cet
  cs_w(r,h)             cost share - welfare
  cs_cl(r,h,*)          cost share - consumption(c)-leisure(l) bundle
  cs_c(t,r,h,*)         cost share - consumption bundle
  cs_cme(t,r,h,*)       cost share - consumption materials and energy bundle
  cs_cene(t,r,h,*)      cost share - consumption electricity and energy bundle
  cs_cen(t,r,h,scen)    cost share - consumption primary energy bundle
  cs_cm(t,r,h,scm)      cost share - consumption materials bundle
  cs_inv(r,s)           cost share - investment bundle
  cs_gov(t,r,s)         cost share - government bundle
  cs_kd_ex(r,s)         cost share - sector differentiated extant capital

* other
  extant_share(r,s)     share of extant capital in production
  lse_sub               labor supply elasticity - substitution
  lse_inc               labor supply elasticity - income
  se_res_sr(s)          short-run supply elasticity for resource sectors
  se_res_lr(s)          long-run supply elasticity for resource sectors
  rate_sr_lr(s)         rate of convergence to long-run supply elasticity for resource sectors
  gov_extra(t,r)        additional government expenditures for use in policy file
  incadj_share(r,h)     household share of additional lump-sum transfers

* stock data
  k0(r)                 benchmark - capital stock
  ke0(r,h)              benchmark - household capital stock
  k_ex0(r)              benchmark - extant capital

* baseline
  n(t,r,h)              baseline - population
  prod_ind(t,r,*,s,*)   baseline - productivity index
  pref_ind(t,r,s,h)     baseline - preference index
  interval              baseline - years between model periods
  beta(t,r,h)           baseline - discount factor
  te(t,r,h)             baseline - time endowment
  cd_ene_growth(s)      baseline - final consumption energy intesity growth rates
  ene_growth(ss,s)      baseline - production energy intensity growth rates
  col_ele_growth        baseline - coal share of electricity production growth rate

* rates
  rho(r,h)              rate - pure rate of time preference
  delta_ex              rate - depreciation rate for extant capital
  delta                 rate - depreciation rate for new capital
  gamma                 rate - population growth rate
  omega                 rate - labor productivity growth rate
  rbar                  rate - benchmark interest rate;



* ------------------------------------------------------------------------------
* define model variables
* ------------------------------------------------------------------------------

positive variables

* prices
  py(t,r,s)             price - output
  pres(t,r,s)           price - fixed factor resource
  pklem(t,r,s)          price - capital(k)-labor(l)-energy(e)-materials(m) bundle
  pkle(t,r,s)           price - capital(k)-labor(l)-energy(e) bundle
  pkl(t,r,s)            price - capital(k)-labor(l) bundle
  pmat(t,r,s)           price - materials(m) bundle
  pene(t,r,s)           price - primary energy(en)-electricity(e) bundle
  pen(t,r,s)            price - primary energy(en) bundle
  pa(t,r,s)             price - armington aggregate
  pdn(t,r,s)            price - domestic(d)-national(n) bundle
  pd(t,r,s)             price - domestic supply
  pn(t,s)               price - national market
  pcl(t,r,h)            price - consumption(c)-leisure(l) bundle
  pc(t,r,h)             price - consumption bundle
  pcme(t,r,h)           price - consumption materials and energy bundle
  pcene(t,r,h)          price - consumption electricity and energy bundle
  pcen(t,r,h)           price - consumption primary energy bundle
  pcm(t,r,h)            price - consumption materials bundle
  pgov(t,r)             price - government bundle
  pfx(t)                price - foreign exchange
  pm(t,s)               price - foreign imports
  px(t,s)               price - foreign exports
  pl(t,r)               price - labor
  pk(t,r)               price - new capital
  pr(t,r)               price - rental rate for new capital
  pr_ex(t,r,s)          price - rental rate for sector specific extant capital
  pr_ex_agg(t,r)        price - rental rate for aggregate extant capital
  pkt(r)                price - terminal capital
  cpi(t)                price - consumer price index

* shadow prices
  lambda(t,r,h)         shadow price - budget constraint

* activities
  y_ex(t,r,s)           activity - output with extant capital
  y(t,r,s)              activity - output with new capital
  a(t,r,s)              activity - armington aggregate
  inv(t,r)              activity - investment bundle
  gov(t,r)              activity - government bundle
  l(t,r,h)              activity - labor

* demands
  ld_ex(t,r,s)          demand - labor for use with extant capital
  kd_ex(t,r,s)          demand - sector specific extant capital
  id_ex(t,r,ss,s)       demand - intermediates for use with extant capital
  res_ex(t,r,s)         demand - fixed factor resource for use with extant capital
  klem(t,r,s)           demand - capital(k)-labor(l)-energy(e)-materials(m) bundle
  kle(t,r,s)            demand - capital(k)-labor(l)-energy(e) bundle
  kl(t,r,s)             demand - capital(k)-labor(l) bundle
  ene(t,r,s)            demand - primary energy(en)-electricity(e) bundle
  en(t,r,s)             demand - primary energy(en) bundle
  mat(t,r,s)            demand - materials(m) bundle
  ld(t,r,s)             demand - labor
  kd(t,r,s)             demand - capital
  id(t,r,ss,s)          demand - intermediates
  res(t,r,s)            demand - fixed factor resource
  d(t,r,s)              demand - domestic supply
  dn(t,r,s)             demand - domestic(d)-national(n) bundle
  m(t,r,s,trd)          demand - imports
  x(t,r,s,trd)          demand - exports
  cl(t,r,h)             demand - consumption(c)-leisure(l) bundle
  leis(t,r,h)           demand - leisure
  c(t,r,h)              demand - consumption bundle
  cme(t,r,h)            demand - consumption materials and energy bundle
  cene(t,r,h)           demand - consumption electricity and energy bundle
  cen(t,r,h)            demand - consumption primary energy bundle
  cm(t,r,h)             demand - consumption materials bundle
  cd(t,r,s,h)           demand - consumption
  i(t,r,s)              demand - investment
  g(t,r,s)              demand - government

* stocks
  k_ex(t,r)             stock - extant capital
  k(t,r)                stock - aggregate new capital
  kh(t,r,h)             stock - household new capital stock
  kt(r)                 stock - terminal capital

* unit costs
  uc_y_ex(t,r,s)        unit cost - output with extant capital
  uc_y(t,r,s)           unit cost - output with new capital
  uc_klem(t,r,s)        unit cost - capital(k)-labor(l)-energy(e)-materials(m) bundle
  uc_kle(t,r,s)         unit cost - capital(k)-labor(l)-energy(e) bundle
  uc_kl(t,r,s)          unit cost - capital(k)-labor(l) bundle
  uc_ene(t,r,s)         unit cost - primary energy(en)-electricity(e) bundle
  uc_en(t,r,s)          unit cost - primary energy(en) bundle
  uc_mat(t,r,s)         unit cost - materials(m) bundle
  uc_a(t,r,s)           unit cost - armington aggregate
  uc_dn(t,r,s)          unit cost - domestic(d)-national(n) armington aggregate
  uc_cl(t,r,h)          unit cost - consumption(c)-leisure(l) bundle
  uc_c(t,r,h)           unit cost - consumption bundle
  uc_cme(t,r,h)         unit cost - consumption materials and energy bundle
  uc_cene(t,r,h)        unit cost - consumption electricity and energy bundle
  uc_cen(t,r,h)         unit cost - consumption primary energy bundle
  uc_cm(t,r,h)          unit cost - consumption materials bundle
  uc_inv(t,r)           unit cost - investment bundle
  uc_gov(t,r)           unit cost - government bundle;

variables

* auxilary
  incadj(t)             auxilary - government to household lump sum transfer
  w(r,h)                auxilary - welfare;



* ------------------------------------------------------------------------------
* define model equations
* ------------------------------------------------------------------------------

equations

* zero profit conditions
  zp_y_ex(t,r,s)        zero profit - output with extant capital
  zp_y(t,r,s)           zero profit - output with new capital
  zp_klem(t,r,s)        zero profit - capital(k)-labor(l)-energy(e)-materials(m) bundle
  zp_kle(t,r,s)         zero profit - capital(k)-labor(l)-energy(e) bundle
  zp_kl(t,r,s)          zero profit - capital(k)-labor(l) bundle
  zp_ene(t,r,s)         zero profit - primary energy(en)-electricity(e) bundle
  zp_en(t,r,s)          zero profit - primary energy(en) bundle
  zp_mat(t,r,s)         zero profit - materials(m) bundle
  zp_a(t,r,s)           zero profit - armington aggregate
  zp_dn(t,r,s)          zero profit - domestic(d)-national(n) armington aggregate
  zp_cl(t,r,h)          zero profit - consumption(c)-leisure(l) bundle
  zp_c(t,r,h)           zero profit - consumption bundle
  zp_cme(t,r,h)         zero profit - consumption materials and energy bundle
  zp_cene(t,r,h)        zero profit - consumption electricity and energy bundle
  zp_cen(t,r,h)         zero profit - consumption primary energy bundle
  zp_cm(t,r,h)          zero profit - consumption materials bundle
  zp_inv(t,r)           zero profit - investment bundle
  zp_k(t,r)             zero profit - capital stock
  zp_gov(t,r)           zero profit - government bundle

* market clearance conditions
  mc_d(t,r,s)           market clearance - domestic output
  mc_a(t,r,s)           market clearance - goods market
  mc_l(t,r)             market clearance - labor
  mc_time(t,r,h)        market clearance - time
  mc_pr_ex(t,r,s)       market clearance - rental rate for sector differentiated extant capital
  mc_pr(t,r)            market clearance - rental rate for new capital
  mc_pk(t,r)            market clearance - new capital
  mc_fx(t)              market clearnace - foreign exchange
  mc_dtrd(t,s)          market clearance - national trade
  mc_pkt(r)             market clearance - terminal capital
  mc_pres(t,r,s)        market clearance - fixed factor resources

* budget constraints
  bc_hh(t,r,h)          budget constraint - household
  bc_gov(t)             budget constraint - government

* hicksian demand equations
  hde_ld_ex(t,r,s)      hicksian demand equation - labor for use with extant capital
  hde_kd_ex(t,r,s)      hicksian demand equation - sector specific extant capital
  hde_id_ex(t,r,ss,s)   hicksian demand equation - intermediates for use with extant capital
  hde_res_ex(t,r,s)     hicksian demand equation - fixed factor resources with extant capital
  hde_klem(t,r,s)       hicksian demand equation - capital(k)-labor(l)-energy(e)-materials(m) bundle
  hde_kle(t,r,s)        hicksian demand equation - capital(k)-labor(l)-energy(e) bundle
  hde_kl(t,r,s)         hicksian demand equation - capital(k)-labor(l) bundle
  hde_ene(t,r,s)        hicksian demand equation - primary energy(en)-electricity(e) bundle
  hde_en(t,r,s)         hicksian demand equation - primary energy(en) bundle
  hde_mat(t,r,s)        hicksian demand equation - materials(m) bundle
  hde_ld(t,r,s)         hicksian demand equation - labor for use with new capital
  hde_kd(t,r,s)         hicksian demand equation - new capital
  hde_id_ele(t,r,sel,s) hicksian demand equation - electricity intermediates
  hde_id_m(t,r,sm,s)    hicksian demand equation - materials intermediates
  hde_id_en(t,r,sen,s)  hicksian demand equation - primary energy intermediates
  hde_res(t,r,s)        hicksian demand equation - fixed factor resource
  hde_m(t,r,s,trd)      hicksian demand equation - imports
  hde_x(t,r,s,trd)      hicksian demand equation - exports
  hde_d(t,r,s)          hicksian demand equation - domestic supply
  hde_dn(t,r,s)         hicksian demand equation - domestic(d)-national(n) bundle
  hde_leis(t,r,h)       hicksian demand equation - leisure
  hde_c(t,r,h)          hicksian demand equation - consumption bundle
  hde_cme(t,r,h)        hicksian demand equation - consumption materials and energy bundle
  hde_cene(t,r,h)       hicksian demand equation - consumption electricity and energy bundle
  hde_cen(t,r,h)        hicksian demand equation - consumption primary energy bundle
  hde_cm(t,r,h)         hicksian demand equation - consumption materials bundle
  hde_cd_trn(t,r,stn,h) hicksian demand equation - consumption transportation
  hde_cd_en(t,r,scen,h) hicksian demand equation - consumption primary energy
  hde_cd_ele(t,r,sel,h) hicksian demand equation - consumption electricity
  hde_cd_m(t,r,scm,h)   hicksian demand equation - consumption materials
  hde_i(t,r,s)          hicksian demand equation - investment
  hde_g(t,r,s)          hicksian demand equation - government

* household first order conditions
  foc_cl(t,r,h)         first order condition - consumption(c)-leisure(l) bundle
  foc_lambda(t,r,h)     first order condition - savings

* unit cost equations
  uce_y_ex(t,r,s)       unit cost equation - output with extant capital
  uce_y(t,r,s)          unit cost equation - output with new capital
  uce_klem(t,r,s)       unit cost equation - capital(k)-labor(l)-energy(e)-materials(m) bundle
  uce_kle(t,r,s)        unit cost equation - capital(k)-labor(l)-energy(e) bundle
  uce_kl(t,r,s)         unit cost equation - capital(k)-labor(l) bundle
  uce_ene(t,r,s)        unit cost equation - primary energy(en)-electricity(e) bundle
  uce_en(t,r,s)         unit cost equation - primary energy(en) bundle
  uce_mat(t,r,s)        unit cost equation - materials(m) bundle
  uce_a(t,r,s)          unit cost equation - armington aggregate
  uce_dn(t,r,s)         unit cost equation - domestic(d)-national(n) armington aggregate
  uce_dx(t,r,s)         unit cost equation - destination differentiated output cet
  uce_k_ex(t,r)         unit cost equation - sector differentiated extant capital
  uce_cl(t,r,h)         unit cost equation - consumption(c)-leisure(l) bundle
  uce_c(t,r,h)          unit cost equation - consumption bundle
  uce_cme(t,r,h)        unit cost equation - consumption materials and energy bundle
  uce_cene(t,r,h)       unit cost equation - consumption electricity and energy bundle
  uce_cen(t,r,h)        unit cost equation - consumption primary energy bundle
  uce_cm(t,r,h)         unit cost equation - consumption materials bundle
  uce_inv(t,r)          unit cost equation - investment bundle
  uce_gov(t,r)          unit cost equation - government bundle

* closure rules
  cr_gov(t,r)           closure rule - real government expenditures fixed
  cr_kt(r)              closure rule - terminal capital

* auxilary
  aux_k_ex(t,r)         auxilary - extant capital
  aux_w(r,h)            auxilary - welfare
  aux_cpi(t)            auxilary - consumer price index;



* ------------------------------------------------------------------------------
* load sets and define subsets
* ------------------------------------------------------------------------------

$gdxin %benchmark_file%
$load r s h trd f

* define subsets
sets
  ftrd(trd)             foreign trade market             / ftrd /
  dtrd(trd)             national trade market            / dtrd /
  sm(s)                 materials sectors                / agf, cru, min, wsu,
                                                           con, fbm, wpm, chm,
                                                           prm, cem, pmm,
                                                           fmm, cpu, tem, bom,
                                                           trn, ttn, srv, hlt /
  sen(s)                primary energy sectors           / col, gas, ref /
  sel(s)                electricity sector               / ele /
  sstd(s)               standard production sectors      / ele, wsu, con, fbm,
                                                           wpm, ref, chm, prm,
                                                           cem, pmm, fmm,
                                                           cpu, tem, bom, trn,
                                                           ttn, srv, hlt /
  sres(s)               resource extraction sectors      / cru, gas, col, min /
  sagf(s)               agriculture and forestry sectors / agf /
  scm(s)                consumption materials sectors    / agf, cru, min, wsu,
                                                           con, fbm, wpm, chm,
                                                           prm, cem, pmm,
                                                           fmm, cpu, tem, bom,
                                                           ttn, srv, hlt /
  stn(s)                transportation sector            / trn /
  scen(s)               consumption energy sectors       / col, gas, ref /;



* ------------------------------------------------------------------------------
* load benchmark data
* ------------------------------------------------------------------------------

* Load benchmark IMPLAN data
$load id0 ld0 kd0 cd0 x0 m0 i0 g0 re0 le0 invh0 n0 ty0 tl0 tc0 tk0



* ------------------------------------------------------------------------------
* load the parameters
* ------------------------------------------------------------------------------

* default optional parameters to zero (i.e., default is to not add features)
res_share(r,s)   = 0;
ene_growth(ss,s) = 0;
cd_ene_growth(s) = 0;
col_ele_growth   = 0;
res_share(r,s)   = 0;
se_res_sr(s)     = 0;
se_res_lr(s)     = 0;
rate_sr_lr(s)    = 0;

* load in parameters
$include %parameter_file%



* ------------------------------------------------------------------------------
* define special temporal sets
* ------------------------------------------------------------------------------

* special sets defining key time periods
tfirst(t) = yes$(ord(t) eq 1);
tlast(t)  = yes$(ord(t) eq card(t));



* ------------------------------------------------------------------------------
* updae rates for the model interval
* ------------------------------------------------------------------------------

* get the interval between model periods
if (card(t) eq 1,
  interval = 1;
else
  interval = sum(t$tfirst(t), sum(tt$(ord(tt) eq (ord(t)+1)),tt.val)-t.val);
);

* convert rates from annual to model time step
gamma              = (1+gamma)**interval-1;
delta              = 1-(1-delta)**interval;
rbar               = (1+rbar)**interval-1;
omega              = (1+omega)**interval-1;

* discount factor via ramsey formula
rho(r,h)           = (1+rbar)/((1+omega+gamma)/(1+gamma))**eta(r,h)-1;
beta(t,r,h)        = (1/(1+rho(r,h)))**(ord(t)-1);



* ------------------------------------------------------------------------------
* disaggregate capital by vintage (extant vs. new capital)
* ------------------------------------------------------------------------------

* if using the putty-clay capital framework disaggregate the benchmark
* activities assuming that the share of capital which is new is based on the
* steady state level of investment - this is only used if there is more than
* one time period in the model
if (%putty_clay% eq 1 and card(t) gt 1,
  extant_share(r,s) = 1-sum(ss, i0(r,ss))/(1+gamma+omega)
                        /(sum(ss, kd0(r,ss))/(rbar+delta));
else
  extant_share(r,s) = 0;);

* don't include the putty-clay specification for resource extraction sectors or
* agriculture and forestry sectors - the fixed factor in those sectors is used
* to calibrate the supply elasticity
extant_share(r,s)$sres(s) = 0;
extant_share(r,s)$sagf(s) = 0;

* set the depreciation rate for extant capital equal to the rate for new capital
delta_ex = delta;

* calculate the inputs into production with extant capital
kd_ex0(r,s)    = extant_share(r,s)*kd0(r,s);
ld_ex0(r,s)    = extant_share(r,s)*ld0(r,s);
id_ex0(r,ss,s) = extant_share(r,s)*id0(r,ss,s);

* calculate the inputs into production with new capital
kd0(r,s)    = (1-extant_share(r,s))*kd0(r,s);
ld0(r,s)    = (1-extant_share(r,s))*ld0(r,s);
id0(r,ss,s) = (1-extant_share(r,s))*id0(r,ss,s);

* production tax is the same whether new or extant capital is used
ty_ex0(r,s) = ty0(r,s);



* ------------------------------------------------------------------------------
* disaggregate land and natural resource payments
* ------------------------------------------------------------------------------

* resource payments for extant capital production
res_ex0(r,s)    = res_share(r,s)*kd_ex0(r,s);
kd_ex0(r,s)     = kd_ex0(r,s)-res_ex0(r,s);
rese_ex0(r,s,h) = res_ex0(r,s)*re0(r,h)/sum(hh, re0(r,hh));

* resource payments for new capital production
res0(r,s)    = res_share(r,s)*kd0(r,s);
kd0(r,s)     = kd0(r,s)-res0(r,s);
rese0(r,s,h) = res0(r,s)*re0(r,h)/sum(hh, re0(r,hh));

* recalculate household capital rental income
re_ex0(r,h) = sum(s, kd_ex0(r,s))*re0(r,h)/sum(hh, re0(r,hh));
re0(r,h)    = sum(s, kd0(r,s))*re0(r,h)/sum(hh, re0(r,hh));



* ------------------------------------------------------------------------------
* baseline cost share adjustments based on energy efficiency improvements
* ------------------------------------------------------------------------------

* initialize the benchmark inputs and final demand
ld_base(t,r,s) = ld0(r,s);
kd_base(t,r,s) = kd0(r,s);
id_base(t,r,ss,s) = id0(r,ss,s);
cd_base(t,r,s,h) = cd0(r,s,h);
g_base(t,r,s) = g0(r,s);

* apply changes in the consumption shares where all changes are balanced by
* sharing out the change across the remaining sectors that were not changed
cd_base(t,r,s,h)$(not cd_ene_growth(s)) = cd_base(t,r,s,h)
                                          +sum(ss$cd_ene_growth(ss),
                                               (1-exp(cd_ene_growth(ss)*(ord(t)-1)*interval))
                                               *cd_base(t,r,ss,h))
                                           *cd_base(t,r,s,h)
                                            /sum(ss$(not cd_ene_growth(ss)), cd_base(t,r,ss,h));
cd_base(t,r,s,h)$cd_ene_growth(s) = exp(cd_ene_growth(s)*(ord(t)-1)*interval)
                                    *cd_base(t,r,s,h);

* for government consumption apply the same energy intensity assumptions as in
* household consumption
g_base(t,r,s)$(not cd_ene_growth(s)) = g_base(t,r,s)
                                      +sum(ss$cd_ene_growth(ss),
                                            (1-exp(cd_ene_growth(ss)*(ord(t)-1)*interval))
                                            *g_base(t,r,ss))
                                       *g_base(t,r,s)
                                        /sum(ss$(not cd_ene_growth(ss)), g_base(t,r,ss));
g_base(t,r,s)$cd_ene_growth(s) = exp(cd_ene_growth(s)*(ord(t)-1)*interval)
                                 *g_base(t,r,s);

* apply energy intesity improvemnts as a change in the cost shares associated
* with energy inputs assuming that the changes are capital embodied resulting
* in a larger capital share
kd_base(t,r,s) = kd_base(t,r,s)
                 +sum(ss, (1-exp(ene_growth(ss,s)*(ord(t)-1)*interval))*id_base(t,r,ss,s))
                  /(1+tk0(r));
id_base(t,r,ss,s) = exp(ene_growth(ss,s)*(ord(t)-1)*interval)*id_base(t,r,ss,s);

* adjust the coal fired electricity generation share assuming that capital and
* labor take up the slack proportionally as a proxy for renewables
kd_base(t,r,s)$(sel(s)) = kd_base(t,r,s)
                          +(1-exp(col_ele_growth*(ord(t)-1)*interval))
                           *id_base(t,r,"col",s)*(1+tk0(r))*kd_base(t,r,s)
                           /((1+tk0(r))*kd_base(t,r,s)+(1+tl0(r))*ld_base(t,r,s))
                           /(1+tk0(r));
ld_base(t,r,s)$(sel(s)) = ld_base(t,r,s)
                          +(1-exp(col_ele_growth*(ord(t)-1)*interval))
                          *id_base(t,r,"col",s)*(1+tl0(r))*ld_base(t,r,s)
                          /((1+tk0(r))*kd_base(t,r,s)+(1+tl0(r))*ld_base(t,r,s))
                          /(1+tl0(r));
id_base(t,r,"col",sel) = exp(col_ele_growth*(ord(t)-1)*interval)*id_base(t,r,"col",sel);



* ------------------------------------------------------------------------------
* calculate remaining benchmark variables
* ------------------------------------------------------------------------------

* commodity bundles in production functions
kl0(t,r,s)           = (1+tl0(r))*ld_base(t,r,s)+(1+tk0(r))*kd_base(t,r,s);
en0(t,r,s)           = sum(sen, id_base(t,r,sen,s));
ene0(t,r,s)          = sum(sel, id_base(t,r,sel,s))+en0(t,r,s);
kle0(t,r,s)          = kl0(t,r,s)+ene0(t,r,s);
mat0(t,r,s)          = sum(sm, id_base(t,r,sm,s));
klem0(t,r,s)         = sum(ss, id_base(t,r,ss,s))+kl0(t,r,s);

* output
y_ex0(r,s)         = ( (1+tl0(r))*ld_ex0(r,s)+(1+tk0(r))*kd_ex0(r,s)
                      +sum(ss, id_ex0(r,ss,s))+(1+tk0(r))*res_ex0(r,s)
                     )/(1-ty_ex0(r,s));
y0(r,s)            = ((1+tl0(r))*ld0(r,s)+(1+tk0(r))*kd0(r,s)
                      +sum(ss, id0(r,ss,s))+(1+tk0(r))*res0(r,s))
                     /(1-ty0(r,s));

* domestic own-use and armington aggregates
d0(r,s)            = y_ex0(r,s)+y0(r,s)-sum(trd, x0(r,s,trd));
dn0(r,s)           = d0(r,s)+m0(r,s,"dtrd");
a0(r,s)            = dn0(r,s)+m0(r,s,"ftrd");

* commodity bundles for utility functions
cen0(t,r,h)          = sum(scen, cd_base(t,r,scen,h))*(1+tc0(r));
cene0(t,r,h)         = cen0(t,r,h)+sum(sel, cd_base(t,r,sel,h))*(1+tc0(r));
cm0(t,r,h)           = sum(scm, cd_base(t,r,scm,h))*(1+tc0(r));
cme0(t,r,h)          = cm0(t,r,h)+cene0(t,r,h);
c0(r,h)              = sum(s, cd0(r,s,h))*(1+tc0(r));

* aggregate regional capitral rents and investment
r0(r)              = sum(h, re0(r,h));
inv0(r)            = sum(s, i0(r,s));

* aggregate regional government expenditures
gov0(r)            = sum(s, g0(r,s));

* balance of payment by household - assumed proportional to consumption
bopdef(r,h)        = (sum(s, m0(r,s,"ftrd")-x0(r,s,"ftrd")))
                     *c0(r,h)/sum(hh, c0(r,hh));
*bopdef(r,h)        = sum((rr,s), m0(rr,s,"ftrd")-x0(rr,s,"ftrd"))
*                     *c0(r,h)/sum((rr,hh), c0(rr,hh));

* benchmark capital stock
k0(r)              = r0(r)/(delta+rbar);
k_ex0(r)           = sum(s, kd_ex0(r,s));

* benchmark capital endowment
ke0(r,h)           = re0(r,h)/(delta+rbar);

* implicit lump sum transfer payments
incadj0(r,h)       = c0(r,h)+invh0(r,h)-re_ex0(r,h)-re0(r,h)-le0(r,h)
                     -sum(s, rese_ex0(r,s,h)+rese0(r,s,h))-bopdef(r,h);

* household labor income
l0(r,h)            = le0(r,h);



* ------------------------------------------------------------------------------
* set consumption-leisure substitution elasticity and benchmark leisure
* ------------------------------------------------------------------------------

* compute benchmark leisure and the c-l substitution elasticity
leis0(r,h)         = -c0(r,h)*lse_inc/(1+lse_inc);
le0(r,h)           = l0(r,h)+leis0(r,h);
se_cl(r,h)         = lse_sub/(1+lse_inc)*(le0(r,h)/leis0(r,h)-1);
cl0(r,h)           = c0(r,h)+leis0(r,h);



* ------------------------------------------------------------------------------
* define additional parameters
* ------------------------------------------------------------------------------

* use benchmark tax rates as defaults
tl(t,r)            = tl0(r);
tk(t,r)            = tk0(r);
tc(t,r)            = tc0(r);
ty_ex(t,r,s)       = ty_ex0(r,s);
ty(t,r,s)          = ty0(r,s);

* household share of additional lump-sum transfers
* lump-sum transfers (incadj) are a national variable used to balance the
* government budget constraint - this parameter determines how that transfer is
* shared out across regions and households
incadj_share(r,h) = c0(r,h)/sum((rr,hh), c0(rr,hh));

* set additional real government spending to zero by default
gov_extra(t,r) = 0;



* ------------------------------------------------------------------------------
* calculate cost shares
* ------------------------------------------------------------------------------

* compute cost shares
cs_y_ex(r,"l",s)        = ((1+tl0(r))*ld_ex0(r,s)          / (y_ex0(r,s)*(1-ty_ex0(r,s))) )$(y_ex0(r,s));
cs_y_ex(r,"k",s)        = ((1+tk0(r))*kd_ex0(r,s)          / (y_ex0(r,s)*(1-ty_ex0(r,s))) )$(y_ex0(r,s));
cs_y_ex(r,ss,s)         = (id_ex0(r,ss,s)                  / (y_ex0(r,s)*(1-ty_ex0(r,s))) )$(y_ex0(r,s));
cs_y_ex(r,"res",s)      = ((1+tk0(r))*res_ex0(r,s)         / (y_ex0(r,s)*(1-ty_ex0(r,s))) )$(y_ex0(r,s));
cs_kd_ex(r,s)           = (kd_ex0(r,s)                     / k_ex0(r)                     )$(k_ex0(r));
cs_rklem(t,r,"res",s)   = ((1+tk0(r))*res0(r,s)            / (y0(r,s)*(1-ty0(r,s)))       )$(y0(r,s));
cs_rklem(t,r,"klem",s)  = (klem0(t,r,s)                    / (y0(r,s)*(1-ty0(r,s)))       )$(y0(r,s));
cs_klem(t,r,"mat",s)    = (mat0(t,r,s)                     / klem0(t,r,s)                 )$(klem0(t,r,s));
cs_klem(t,r,"kle",s)    = (kle0(t,r,s)                     / klem0(t,r,s)                 )$(klem0(t,r,s));
cs_kle(t,r,"ene",s)     = (ene0(t,r,s)                     / kle0(t,r,s)                  )$(kle0(t,r,s));
cs_kle(t,r,"kl",s)      = (kl0(t,r,s)                      / kle0(t,r,s)                  )$(kle0(t,r,s));
cs_kl(t,r,"k",s)        = ((1+tk0(r))*kd_base(t,r,s)       / kl0(t,r,s)                   )$(kl0(t,r,s));
cs_kl(t,r,"l",s)        = ((1+tl0(r))*ld_base(t,r,s)       / kl0(t,r,s)                   )$(kl0(t,r,s));
cs_ene(t,r,sel,s)       = (id_base(t,r,sel,s)              / ene0(t,r,s)                  )$(ene0(t,r,s));
cs_ene(t,r,"en",s)      = (en0(t,r,s)                      / ene0(t,r,s)                  )$(ene0(t,r,s));
cs_en(t,r,sen,s)        = (id_base(t,r,sen,s)              / en0(t,r,s)                   )$(en0(t,r,s));
cs_mat(t,r,sm,s)        = (id_base(t,r,sm,s)               / mat0(t,r,s)                  )$(mat0(t,r,s));
cs_dn(r,"d",s)          = (d0(r,s)                         / dn0(r,s)                     )$(dn0(r,s));
cs_dn(r,"n",s)          = (m0(r,s,"dtrd")                  / dn0(r,s)                     )$(dn0(r,s));
cs_nf(r,"n",s)          = (dn0(r,s)                        / a0(r,s)                      )$(a0(r,s));
cs_nf(r,"f",s)          = (m0(r,s,"ftrd")                  / a0(r,s)                      )$(a0(r,s));
cs_dx(r,"d",s)          = (d0(r,s)                         / (y_ex0(r,s)+y0(r,s))         )$(y_ex0(r,s)+y0(r,s));
cs_dx(r,"n",s)          = (x0(r,s,"dtrd")                  / (y_ex0(r,s)+y0(r,s))         )$(y_ex0(r,s)+y0(r,s));
cs_dx(r,"f",s)          = (x0(r,s,"ftrd")                  / (y_ex0(r,s)+y0(r,s))         )$(y_ex0(r,s)+y0(r,s));
cs_c(t,r,h,"trn")       = ((1+tc0(r))*cd_base(t,r,"trn",h) / c0(r,h)                      )$(c0(r,h));
cs_c(t,r,h,"cme")       = (cme0(t,r,h)                     / c0(r,h)                      )$(c0(r,h));
cs_cme(t,r,h,"cm")      = (cm0(t,r,h)                      / cme0(t,r,h)                  )$(cme0(t,r,h));
cs_cme(t,r,h,"cene")    = (cene0(t,r,h)                    / cme0(t,r,h)                  )$(cme0(t,r,h));
cs_cene(t,r,h,"cen")    = (cen0(t,r,h)                     / cene0(t,r,h)                 )$(cene0(t,r,h));
cs_cene(t,r,h,"ele")    = ((1+tc0(r))*cd_base(t,r,"ele",h) / cene0(t,r,h)                 )$(cene0(t,r,h));
cs_cen(t,r,h,scen)      = ((1+tc0(r))*cd_base(t,r,scen,h)  / cen0(t,r,h)                  )$(cen0(t,r,h));
cs_cm(t,r,h,scm)        = ((1+tc0(r))*cd_base(t,r,scm,h)   / cm0(t,r,h)                   )$(cm0(t,r,h));
cs_cl(r,h,"leis")       = (leis0(r,h)                      / cl0(r,h)                     )$(cl0(r,h));
cs_cl(r,h,"c")          = (c0(r,h)                         / cl0(r,h)                     )$(cl0(r,h));
cs_inv(r,s)             = (i0(r,s)                         / inv0(r)                      )$(inv0(r));
cs_gov(t,r,s)           = (g_base(t,r,s)                   / gov0(r)                      )$(gov0(r));



* ------------------------------------------------------------------------------
* calibrate resource substitution elasticity in extraction sectors
* ------------------------------------------------------------------------------

* substitution elasticity - r-klem bundle in resource extraction sectors
* for the dynamic model the substitution elasticity is calibrated such that the
* initial price elasticity of supply for the sector is equal to se_res_sr
* but where the elasticity of supply convergences at rate rate_sr_lr to the
* price elasticity of supply se_res_lr
* for the static model we calibrate the price elasticity of supply to be equal
* to se_res_lr
if (card(t) eq 1,
  se_rklem(t,r,s)$cs_rklem(t,r,"res",s) = se_res_lr(s)/(1/cs_rklem(t,r,"res",s)-1);
else
  se_rklem(t,r,s)$cs_rklem(t,r,"res",s) = (se_res_lr(s)-(se_res_lr(s)-se_res_sr(s))
                                                        *exp(-rate_sr_lr(s)
                                                              *(ord(t)-1)*interval)
                                          )/(1/cs_rklem(t,r,"res",s)-1);
);


* ------------------------------------------------------------------------------
* define the baseline
* ------------------------------------------------------------------------------

* baseline household population grows at rate gamma
n(t,r,h)  = (1+gamma)**(ord(t)-1)*n0(r,h);

* time endowment based on population growth and harrod neutral growth rate
te(t,r,h) = (1+gamma+omega)**(ord(t)-1)*le0(r,h);

* initialize productivity index to unity
prod_ind(t,r,ss,s,"new")      = 1;
prod_ind(t,r,ss,s,"extant")   = 1;
prod_ind(t,r,"k",s,"new")     = 1;
prod_ind(t,r,"k",s,"extant")  = 1;
prod_ind(t,r,"l",s,"new")     = 1;
prod_ind(t,r,"l",s,"extant")  = 1;

* intialize preference index to unity
pref_ind(t,r,s,h) = 1;



* ------------------------------------------------------------------------------
* unit cost functions
* ------------------------------------------------------------------------------

* unit cost equation - output with extant capital
uce_y_ex(t,r,s)$(y_ex0(r,s))..
  uc_y_ex(t,r,s) =e= (sum(ss, cs_y_ex(r,ss,s)*prod_ind(t,r,ss,s,"extant")*pa(t,r,ss))
                      +cs_y_ex(r,"l",s)*prod_ind(t,r,"l",s,"extant")*(1+tl(t,r))*pl(t,r)/(1+tl0(r))
                      +cs_y_ex(r,"k",s)*prod_ind(t,r,"k",s,"extant")*(1+tk(t,r))*pr_ex(t,r,s)/(1+tk0(r))
                      +cs_y_ex(r,"res",s)*(1+tk(t,r))*pres(t,r,s)/(1+tk0(r))
                     )*(1-ty_ex0(r,s));

* unit cost equation - output with new capital
uce_y(t,r,s)$(y0(r,s))..
  uc_y(t,r,s) =e= ( cs_rklem(t,r,"res", s)*((1+tk(t,r))*pres(t,r,s)/(1+tk0(r)))**(1-se_rklem(t,r,s))
                   +cs_rklem(t,r,"klem",s)*pklem(t,r,s)                        **(1-se_rklem(t,r,s))
                  )**(1/(1-se_rklem(t,r,s)))*(1-ty0(r,s));

* unit cost equation - capital(k)-labor(l)-energy(e)-materials(m) bundle
uce_klem(t,r,s)$(y0(r,s))..
  uc_klem(t,r,s) =e= ( cs_klem(t,r,"mat",s)*pmat(t,r,s)**(1-se_klem(s))
                      +cs_klem(t,r,"kle",s)*pkle(t,r,s)**(1-se_klem(s))
                     )**(1/(1-se_klem(s)));

* unit cost equation - capital(k)-labor(l)-energy(e) bundle
uce_kle(t,r,s)$(y0(r,s))..
  uc_kle(t,r,s) =e= ( cs_kle(t,r,"ene",s)*pene(t,r,s)**(1-se_kle(s))
                     +cs_kle(t,r,"kl",s) *pkl(t,r,s) **(1-se_kle(s))
                    )**(1/(1-se_kle(s)));

* unit cost equation - capital(k)-labor(l) bundle
uce_kl(t,r,s)$(y0(r,s))..
  uc_kl(t,r,s) =e= ( cs_kl(t,r,"l",s)*(prod_ind(t,r,"l",s,"new")*(1+tl(t,r))*pl(t,r)/(1+tl0(r)))**(1-se_kl(s))
                    +cs_kl(t,r,"k",s)*(prod_ind(t,r,"k",s,"new")*(1+tk(t,r))*pr(t,r)/(1+tk0(r)))**(1-se_kl(s))
                   )**(1/(1-se_kl(s)));

* unit cost equation - primary energy(en)-electricity(e) bundle
uce_ene(t,r,s)$(y0(r,s))..
  uc_ene(t,r,s) =e= ( cs_ene(t,r,"ele",s)*(sum(ss$sel(ss), prod_ind(t,r,ss,s,"new")*pa(t,r,ss)))**(1-se_ene(s))
                     +cs_ene(t,r,"en",s) *(pen(t,r,s))                                          **(1-se_ene(s))
                    )**(1/(1-se_ene(s)));

* unit cost equation - primary energy(en) bundle
uce_en(t,r,s)$(y0(r,s))..
  uc_en(t,r,s) =e= sum(sen, cs_en(t,r,sen,s)*(prod_ind(t,r,sen,s,"new")*pa(t,r,sen))**(1-se_en(s)))
                   **(1/(1-se_en(s)));

* unit cost equation - materials(m) bundle
uce_mat(t,r,s)$(y0(r,s))..
  uc_mat(t,r,s) =e= sum(sm, cs_mat(t,r,sm,s)*prod_ind(t,r,sm,s,"new")*pa(t,r,sm));

* unit cost equation - armington aggregate
uce_a(t,r,s)..
  uc_a(t,r,s) =e= ( cs_nf(r,"n",s)*pdn(t,r,s)      **(1-se_nf(s))
                   +cs_nf(r,"f",s)*(pfx(t)*pm(t,s))**(1-se_nf(s))
                  )**(1/(1-se_nf(s)));

* unit cost equation - domestic-national armington aggregate
uce_dn(t,r,s)..
  uc_dn(t,r,s) =e= ( (cs_dn(r,"d",s)*pd(t,r,s)**(1-se_dn(s)))$(cs_dn(r,"d",s))
                    +(cs_dn(r,"n",s)*pn(t,s)  **(1-se_dn(s)))$(cs_dn(r,"n",s))
                   )**(1/(1-se_dn(s)));

* unit cost equation - destination differentiated output cet
uce_dx(t,r,s)$(y0(r,s))..
  py(t,r,s) =e= ( cs_dx(r,"d",s)*pd(t,r,s)       **(1+te_dx(s))
                 +cs_dx(r,"n",s)*pn(t,s)         **(1+te_dx(s))
                 +cs_dx(r,"f",s)*(pfx(t)*px(t,s))**(1+te_dx(s))
                )**(1/(1+te_dx(s)));

* unit cost equation - sector differentiated extant capital
uce_k_ex(t,r)$(k_ex0(r))..
  pr_ex_agg(t,r) =e= sum(s, cs_kd_ex(r,s)*pr_ex(t,r,s)**(1+te_k_ex))
                     **(1/(1+te_k_ex));

* unit cost equation - consumption(c)-leisure(l) bundle
uce_cl(t,r,h)..
  uc_cl(t,r,h) =e= ( cs_cl(r,h,"leis")*pl(t,r)  **(1-se_cl(r,h))
                    +cs_cl(r,h,"c"   )*pc(t,r,h)**(1-se_cl(r,h))
                   )**(1/(1-se_cl(r,h)));

* unit cost equation - consumption bundle
uce_c(t,r,h)..
  uc_c(t,r,h) =e= ( cs_c(t,r,h,"trn")*((1+tc(t,r))
                                     *sum(s$stn(s), pref_ind(t,r,"trn",h)
                                                    *pa(t,r,s)/(1+tc0(r)))
                                    )**(1-se_c)
                   +cs_c(t,r,h,"cme")*pcme(t,r,h)**(1-se_c)
                  )**(1/(1-se_c));

* unit cost equation - consumption materials and energy bundle
uce_cme(t,r,h)..
  uc_cme(t,r,h) =e= ( cs_cme(t,r,h,"cm"  )*pcm(t,r,h)  **(1-se_cme)
                     +cs_cme(t,r,h,"cene")*pcene(t,r,h)**(1-se_cme)
                    )**(1/(1-se_cme));

* unit cost equation - consumption electricity and energy bundle
uce_cene(t,r,h)..
  uc_cene(t,r,h) =e= ( cs_cene(t,r,h,"ele")*((1+tc(t,r))
                                           *sum(s$sel(s), pref_ind(t,r,"ele",h)
                                                          *pa(t,r,s)/(1+tc0(r)))
                                          )**(1-se_cene)
                      +cs_cene(t,r,h,"cen")*pcen(t,r,h)**(1-se_cene)
                     )**(1/(1-se_cene));

* unit cost equation - consumption primary energy bundle
uce_cen(t,r,h)..
  uc_cen(t,r,h) =e= sum(scen, cs_cen(t,r,h,scen)*(pref_ind(t,r,scen,h)
                                                *(1+tc(t,r))*pa(t,r,scen)/(1+tc0(r))
                                               )**(1-se_cen))**(1/(1-se_cen));

* unit cost equation - consumption materials bundle
uce_cm(t,r,h)..
  uc_cm(t,r,h) =e= sum(scm, cs_cm(t,r,h,scm)*(pref_ind(t,r,scm,h)
                                            *(1+tc(t,r))*pa(t,r,scm)/(1+tc0(r))
                                           )**(1-se_cm))**(1/(1-se_cm));

* unit cost equation - investment bundle
uce_inv(t,r)..
  uc_inv(t,r) =e= sum(s, cs_inv(r,s)*pa(t,r,s));

* unit cost equation - government bundle
uce_gov(t,r)..
  uc_gov(t,r) =e= sum(s, cs_gov(t,r,s)*pa(t,r,s));



* ------------------------------------------------------------------------------
* hicksian demand functions - via shepard's lemma
* ------------------------------------------------------------------------------

* hicksian demand equation - labor for use with extant capital
hde_ld_ex(t,r,s)$(ld_ex0(r,s))..
  ld_ex(t,r,s) =e= prod_ind(t,r,"l",s,"extant")*y_ex(t,r,s);

* hicksian demand equation - sector specific extant capital
hde_kd_ex(t,r,s)$(kd_ex0(r,s))..
  kd_ex(t,r,s) =e= prod_ind(t,r,"k",s,"extant")*y_ex(t,r,s);

* hicksian demand equation - intermediates for use with extant capital
hde_id_ex(t,r,ss,s)$(id_ex0(r,ss,s))..
  id_ex(t,r,ss,s) =e= prod_ind(t,r,ss,s,"extant")*y_ex(t,r,s);

* hicksian demand equation - fixed factor resources with extant capital
hde_res_ex(t,r,s)$(res_ex0(r,s))..
  res_ex(t,r,s) =e= y_ex(t,r,s);

* hicksian demand equation - fixed factor resource
hde_res(t,r,s)$(res0(r,s))..
  res(t,r,s) =e= y(t,r,s)*( uc_y(t,r,s)/(1-ty0(r,s))
                           *(1+tk0(r))/((1+tk(t,r))*pres(t,r,s))
                          )**se_rklem(t,r,s);

* hicksian demand equation - capital(k)-labor(l)-energy(e)-materials(m) bundle
hde_klem(t,r,s)$(y0(r,s))..
  klem(t,r,s) =e= y(t,r,s)*(uc_y(t,r,s)/(1-ty0(r,s))/pklem(t,r,s))
                           **se_rklem(t,r,s);

* hicksian demand equation - capital(k)-labor(l)-energy(e) bundle
hde_kle(t,r,s)$(y0(r,s))..
  kle(t,r,s) =e= klem(t,r,s)*(uc_klem(t,r,s)/pkle(t,r,s))**se_klem(s);

* hicksian demand equation - capital(k)-labor(l) bundle
hde_kl(t,r,s)$(y0(r,s))..
  kl(t,r,s) =e= kle(t,r,s)*(uc_kle(t,r,s)/pkl(t,r,s))**se_kle(s);

* hicksian demand equation - labor for use with new capital
hde_ld(t,r,s)$(y0(r,s))..
  ld(t,r,s) =e= kl(t,r,s)*prod_ind(t,r,"l",s,"new")**(1-se_kl(s))
                *(uc_kl(t,r,s)*(1+tl0(r))/((1+tl(t,r))*pl(t,r)))**se_kl(s);

* hicksian demand equation - new capital
hde_kd(t,r,s)$(y0(r,s))..
  kd(t,r,s) =e= kl(t,r,s)*prod_ind(t,r,"k",s,"new")**(1-se_kl(s))
                *(uc_kl(t,r,s)*(1+tk0(r))/((1+tk(t,r))*pr(t,r)))**se_kl(s);

* hicksian demand equation - primary energy(en)-electricity(e) bundle
hde_ene(t,r,s)$(y0(r,s))..
  ene(t,r,s) =e= kle(t,r,s)*(uc_kle(t,r,s)/pene(t,r,s))**se_kle(s);

* hicksian demand equation - electricity intermediates
hde_id_ele(t,r,sel,s)$(y0(r,s))..
  id(t,r,sel,s) =e= ene(t,r,s)*prod_ind(t,r,sel,s,"new")**(1-se_ene(s))
                      *(uc_ene(t,r,s)/pa(t,r,sel))**se_ene(s);

* hicksian demand equation - primary energy(en) bundle
hde_en(t,r,s)$(y0(r,s))..
  en(t,r,s) =e= ene(t,r,s)*(uc_ene(t,r,s)/pen(t,r,s))**se_ene(s);

* hicksian demand equation - materials(m) bundle
hde_mat(t,r,s)$(y0(r,s))..
  mat(t,r,s) =e= klem(t,r,s)*(uc_klem(t,r,s)/pmat(t,r,s))**se_klem(s);

* hicksian demand equation - materials intermediates
hde_id_m(t,r,sm,s)$(y0(r,s))..
  id(t,r,sm,s) =e= prod_ind(t,r,sm,s,"new")*mat(t,r,s);

* hicksian demand equation - primary energy
hde_id_en(t,r,sen,s)$(y0(r,s))..
  id(t,r,sen,s) =e= en(t,r,s)*prod_ind(t,r,sen,s,"new")**(1-se_en(s))
                    *(uc_en(t,r,s)/pa(t,r,sen))**se_en(s);

* hicksian demand equation - domestic-national armington aggregate
hde_dn(t,r,s)..
  dn(t,r,s) =e= a(t,r,s)*(uc_a(t,r,s)/pdn(t,r,s))**se_nf(s);

* hicksian demand equation - imports
hde_m(t,r,s,trd)$(m0(r,s,trd))..
  m(t,r,s,trd) =e=  (dn(t,r,s)*(uc_dn(t,r,s)/pn(t,s))        **se_dn(s))$(dtrd(trd))
                   +( a(t,r,s)*(uc_a(t,r,s)/(pfx(t)*pm(t,s)))**se_nf(s))$(ftrd(trd));

* hicksian demand equation - domestic supply
hde_d(t,r,s)$(d0(r,s))..
  d(t,r,s) =e= dn(t,r,s)*(uc_dn(t,r,s)/pd(t,r,s))**se_dn(s);

* hicksian demand equation - exports
hde_x(t,r,s,trd)$(x0(r,s,trd))..
  x(t,r,s,trd) =e= (y_ex(t,r,s)*y_ex0(r,s)+y(t,r,s)*y0(r,s))/(y_ex0(r,s)+y0(r,s))
                   *( ((pn(t,s)/py(t,r,s))       **te_dx(s))$(dtrd(trd))
                     +((pfx(t)*px(t,s)/py(t,r,s))**te_dx(s))$(ftrd(trd)));

* hicksian demand equation - leisure
hde_leis(t,r,h)..
  leis(t,r,h) =e= cl(t,r,h)*(uc_cl(t,r,h)/pl(t,r))**se_cl(r,h);

* hicksian demand equation - consumption bundle
hde_c(t,r,h)..
  c(t,r,h) =e= cl(t,r,h)*(uc_cl(t,r,h)/pc(t,r,h))**se_cl(r,h);

* hicksian demand equation - consumption materials and energy bundle
hde_cme(t,r,h)..
  cme(t,r,h) =e= c(t,r,h)*(uc_c(t,r,h)/pcme(t,r,h))**se_c;

* hicksian demand equation - consumption transportation
hde_cd_trn(t,r,stn,h)..
  cd(t,r,stn,h) =e= c(t,r,h)*pref_ind(t,r,stn,h)**(1-se_c)
                    *(uc_c(t,r,h)*(1+tc0(r))/((1+tc(t,r))*pa(t,r,stn)))**se_c;

* hicksian demand equation - consumption materials bundle
hde_cm(t,r,h)..
  cm(t,r,h) =e= cme(t,r,h)*(uc_cme(t,r,h)/pcm(t,r,h))**se_cme;

* hicksian demand equation - consumption electricity and energy bundle
hde_cene(t,r,h)..
  cene(t,r,h) =e= cme(t,r,h)*(uc_cme(t,r,h)/pcene(t,r,h))**se_cme;

* hicksian demand equation - consumption electricity
hde_cd_ele(t,r,sel,h)..
  cd(t,r,sel,h) =e= cene(t,r,h)*pref_ind(t,r,sel,h)**(1-se_cene)
                    *(uc_cene(t,r,h)*(1+tc0(r))/((1+tc(t,r))*pa(t,r,sel)))**se_cene;

* hicksian demand equation - consumption primary energy bundle
hde_cen(t,r,h)..
  cen(t,r,h) =e= cene(t,r,h)*(uc_cene(t,r,h)/pcen(t,r,h))**se_cene;

* hicksian demand equation - consumption materials
hde_cd_m(t,r,scm,h)..
  cd(t,r,scm,h) =e= cm(t,r,h)*pref_ind(t,r,scm,h)**(1-se_cm)
                    *(uc_cm(t,r,h)*(1+tc0(r))/((1+tc(t,r))*pa(t,r,scm)))**se_cm;

* hicksian demand equation - consumption energy
hde_cd_en(t,r,scen,h)..
  cd(t,r,scen,h) =e= cen(t,r,h)*pref_ind(t,r,scen,h)**(1-se_cen)
                     *(uc_cen(t,r,h)*(1+tc0(r))/((1+tc(t,r))*pa(t,r,scen)))**se_cen;

* hicksian demand equation - investment
hde_i(t,r,s)..
  i(t,r,s) =e= inv(t,r);

* hicksian demand equation - government
hde_g(t,r,s)..
  g(t,r,s) =e= gov(t,r);



* ------------------------------------------------------------------------------
* household first order conditions
* ------------------------------------------------------------------------------

* first order condition - consumption(c)-leisure(l) bundle
foc_cl(t,r,h)..
  (cl(t,r,h)*cl0(r,h)/n(t,r,h))**(-eta(r,h)) =e= lambda(t,r,h)*pcl(t,r,h);

* first order condition - savings
foc_lambda(t+1,r,h)..
  beta(t+1,r,h)*lambda(t+1,r,h) =e= beta(t,r,h)*lambda(t,r,h);



* ------------------------------------------------------------------------------
* zero profit conditions
* ------------------------------------------------------------------------------

* zero profit - output with extant capital
zp_y_ex(t,r,s)$(y_ex0(r,s))..
  uc_y_ex(t,r,s) =e= py(t,r,s)*(1-ty_ex(t,r,s));

* zero profit - output with new capital
zp_y(t,r,s)$(y0(r,s))..
  uc_y(t,r,s) =e= py(t,r,s)*(1-ty(t,r,s));

* zero profit - capital(k)-labor(l)-energy(e)-materials(m) bundle
zp_klem(t,r,s)$(y0(r,s))..
  uc_klem(t,r,s) =e= pklem(t,r,s);

* zero profit - capital(k)-labor(l)-energy(e) bundle
zp_kle(t,r,s)$(y0(r,s))..
  uc_kle(t,r,s) =e= pkle(t,r,s);

* zero profit - capital(k)-labor(l) bundle
zp_kl(t,r,s)$(y0(r,s))..
  uc_kl(t,r,s) =e= pkl(t,r,s);

* zero profit - primary energy(en)-electricity(e) bundle
zp_ene(t,r,s)$(y0(r,s))..
  uc_ene(t,r,s) =e= pene(t,r,s);

* zero profit - primary energy(en) bundle
zp_en(t,r,s)$(y0(r,s))..
  uc_en(t,r,s) =e= pen(t,r,s);

* zero profit - materials(m) bundle
zp_mat(t,r,s)$(y0(r,s))..
  uc_mat(t,r,s) =e= pmat(t,r,s);

* zero profit - armington aggregate
zp_a(t,r,s)$(a0(r,s))..
  uc_a(t,r,s) =e= pa(t,r,s);

* zero profit - domestic-national armington aggregate
zp_dn(t,r,s)$(dn0(r,s))..
  uc_dn(t,r,s) =e= pdn(t,r,s);

* zero profit - consumption-leisure bundle
zp_cl(t,r,h)$(cl0(r,h))..
  uc_cl(t,r,h) =e= pcl(t,r,h);

* zero profit - consumption bundle
zp_c(t,r,h)$(c0(r,h))..
  uc_c(t,r,h) =e= pc(t,r,h);

* zero profit - consumption materials and energy bundle
zp_cme(t,r,h)$(cme0(t,r,h))..
  uc_cme(t,r,h) =e= pcme(t,r,h);

* zero profit - consumption electricity and energy bundle
zp_cene(t,r,h)$(cene0(t,r,h))..
  uc_cene(t,r,h) =e= pcene(t,r,h);

* zero profit - consumption primary energy bundle
zp_cen(t,r,h)$(cen0(t,r,h))..
  uc_cen(t,r,h) =e= pcen(t,r,h);

* zero profit - consumption materials bundle
zp_cm(t,r,h)$(cm0(t,r,h))..
  uc_cm(t,r,h) =e= pcm(t,r,h);

* zero profit - investment
zp_inv(t,r)$(inv0(r))..
  uc_inv(t,r) =e= pk(t+1,r)+pkt(r)$tlast(t);

* zero profit - capital stock
zp_k(t,r)..
  pk(t,r) =e= pr(t,r)*(rbar+delta)+(pk(t+1,r)+pkt(r)$tlast(t))*(1-delta);

* zero profit - government bundle
zp_gov(t,r)$(gov0(r))..
  uc_gov(t,r) =e= pgov(t,r);



* ------------------------------------------------------------------------------
* market clearance conditions
* ------------------------------------------------------------------------------

* market clearance - domestic output
mc_d(t,r,s)$(y_ex0(r,s)+y0(r,s))..
  (y_ex(t,r,s)*y_ex0(r,s)+y(t,r,s)*y0(r,s))/(y_ex0(r,s)+y0(r,s))
    *(pd(t,r,s)/py(t,r,s))**te_dx(s)
      =e= d(t,r,s);

* market clearance - goods market
mc_a(t,r,s)$(a0(r,s))..
  a(t,r,s)*a0(r,s) =e=  sum(ss, id(t,r,s,ss)*id_base(t,r,s,ss))
                       +sum(ss, id_ex(t,r,s,ss)*id_ex0(r,s,ss))
                       +sum(h, cd(t,r,s,h)*cd_base(t,r,s,h))
                       +i(t,r,s)*i0(r,s)+g(t,r,s)*g_base(t,r,s);

* market clearance - labor
mc_l(t,r)..
  sum(h, l(t,r,h)*l0(r,h)) =e=  sum(s, ld_ex(t,r,s)*ld_ex0(r,s))
                               +sum(s, ld(t,r,s)*ld_base(t,r,s));

* market clearance - time
mc_time(t,r,h)..
 te(t,r,h) =e= leis(t,r,h)*leis0(r,h)+l(t,r,h)*l0(r,h);

* market clearance - rental rate for sector differentiated extant capital
mc_pr_ex(t,r,s)$(kd_ex0(r,s))..
  k_ex(t,r)*(pr_ex(t,r,s)/pr_ex_agg(t,r))**te_k_ex =e= kd_ex(t,r,s);

* market clearance - rental rate for new capital
mc_pr(t,r)..
  k(t,r)*r0(r) =e= sum(s, kd(t,r,s)*kd_base(t,r,s));

* market clearance - price of new capital
mc_pk(t,r)..
  k(t-1,r)*k0(r)*(1-delta)+inv(t-1,r)*inv0(r)+k0(r)$tfirst(t)
    =e= k(t,r)*k0(r);

* market clearnace - foreign trade
mc_fx(t)..
  sum((r,s,ftrd), px(t,s)*x(t,r,s,ftrd)*x0(r,s,ftrd))
  +sum((r,h), bopdef(r,h)*(1+gamma+omega)**(ord(t)-1))
    =e= sum((r,s,ftrd), pm(t,s)*m(t,r,s,ftrd)*m0(r,s,ftrd));

* market clearance - national trade
mc_dtrd(t,s)..
  sum((r,dtrd), x(t,r,s,dtrd)*x0(r,s,dtrd))
    =e= sum((r,dtrd), m(t,r,s,dtrd)*m0(r,s,dtrd));

* market clearance - price of post-terminal capital
mc_pkt(r)..
  sum(t$tlast(t), k(t,r)*k0(r)*(1-delta)+inv(t,r)*inv0(r)) =e= kt(r)*k0(r);

* market clearnace - fixed factor resources
mc_pres(t,r,s)..
  res_ex0(r,s)+res0(r,s) =g= res_ex(t,r,s)*res_ex0(r,s)+res(t,r,s)*res0(r,s);



* ------------------------------------------------------------------------------
* budget constraints
* ------------------------------------------------------------------------------

* budget constraint - household
* the handling of the capital endowment using the conditionals based on the
* length of the t set is to correctly specify the budget constraint for the
* static model where the value of the endowment is defined by pr and not pk
bc_hh(t,r,h)..
  pk(t+1,r)*kh(t+1,r,h)*ke0(r,h)
  +pkt(r)*kt(r)*ke0(r,h)$(tlast(t) and card(t) gt 1)
  +pkt(r)*invh0(r,h)$(card(t) eq 1)
  +pcl(t,r,h)*cl(t,r,h)*cl0(r,h)
    =e=  pk(t,r)*kh(t,r,h)*ke0(r,h)$(card(t) gt 1)
         +pr(t,r)*ke0(r,h)$(card(t) eq 1)
         +pl(t,r)*te(t,r,h)
         +pr_ex_agg(t,r)*k_ex(t,r)*re_ex0(r,h)
         +sum(s, pres(t,r,s)*( res_ex(t,r,s)*rese_ex0(r,s,h)
                              +res(t,r,s)*rese0(r,s,h)))
          +pfx(t)*bopdef(r,h)*(1+gamma+omega)**(ord(t)-1)
          +cpi(t)*incadj0(r,h)*(1+gamma+omega)**(ord(t)-1)
          +cpi(t)*incadj(t)*incadj_share(r,h);

* budget constraint - government
bc_gov(t)..
  sum(r, pgov(t,r)*gov(t,r)*gov0(r))
    =e= sum(r, sum(s, tl(t,r)*pl(t,r)*ld(t,r,s)*ld_base(t,r,s)
                     +tl(t,r)*pl(t,r)*ld_ex(t,r,s)*ld_ex0(r,s)
                     +tk(t,r)*pr(t,r)*kd(t,r,s)*kd_base(t,r,s)
                     +tk(t,r)*pres(t,r,s)*( res_ex(t,r,s)*res_ex0(r,s)
                                           +res(t,r,s)*res0(r,s))
                     +tk(t,r)*pr_ex(t,r,s)*kd_ex(t,r,s)*kd_ex0(r,s)
                     +ty_ex(t,r,s)*py(t,r,s)*y_ex(t,r,s)*y_ex0(r,s)
                     +ty(t,r,s)*py(t,r,s)*y(t,r,s)*y0(r,s)
                     +sum(h, tc(t,r)*pa(t,r,s)*cd(t,r,s,h)*cd_base(t,r,s,h)))
                     -sum(h, cpi(t)*incadj0(r,h)*(1+gamma+omega)**(ord(t)-1)))
        -cpi(t)*incadj(t);



* ------------------------------------------------------------------------------
* closure rules
* ------------------------------------------------------------------------------

* closure rule - real government expenditures per capita are fixed
cr_gov(t,r)..
  gov(t,r)*gov0(r) =e= (1+gamma+omega)**(ord(t)-1)*gov0(r)+gov_extra(t,r);
*  gov(t,r)*gov0(r) =e= (1+gamma+omega)**(ord(t)-1)*gov0(r)+gov_extra(t,r);

* closure rule - terminal capital
* scale the level of the terminal capital stock to achieve steady-state growth
* in last period investment.
cr_kt(r)..
  sum(tlast(t), inv(t,r)*sum(h, c(t-1,r,h)))
    =e= sum(tlast(t), sum(h, c(t,r,h))*inv(t-1,r));



* ------------------------------------------------------------------------------
* auxilary
* ------------------------------------------------------------------------------

* auxilary - extant capital
aux_k_ex(t,r)$(k_ex0(r))..
  k_ex(t,r) =e= (1-delta_ex)**(ord(t)-1);

* auxilary - welfare
aux_w(r,h)..
  w(r,h) =e= sum(t,  n(t,r,h)*beta(t,r,h)
                    *(cl(t,r,h)*cl0(r,h)/n(t,r,h))**(1-eta(r,h))/(1-eta(r,h)));

* auxilary - consumer price index
aux_cpi(t)..
  cpi(t) =e= sum((rr,hh), pc(t,rr,hh)*c0(rr,hh))/sum((rr,hh), c0(rr,hh));



* ------------------------------------------------------------------------------
* define mcp problem
* ------------------------------------------------------------------------------

model sage /

* zero profit conditions - determine demand for bundled commodities
  zp_y_ex.y_ex
  zp_y.y
  zp_klem.klem
  zp_kle.kle
  zp_kl.kl
  zp_ene.ene
  zp_en.en
  zp_mat.mat
  zp_a.a
  zp_dn.dn
  zp_cl.pcl
  zp_c.c
  zp_cme.cme
  zp_cene.cene
  zp_cen.cen
  zp_cm.cm
  zp_inv.inv
  zp_k.k
  zp_gov.gov

* market clearance conditions - determine prices for underlying commodities
  mc_pr_ex.pr_ex
  mc_pres.pres
  mc_d.pd
  mc_a.pa
  mc_pr.pr
  mc_pk.pk
  mc_fx.pfx
  mc_pkt.pkt
  mc_time.l
  mc_l.pl
  mc_dtrd.pn

* hicksian demand equations - determine prices for bundled commodities
  hde_klem.pklem
  hde_kle.pkle
  hde_kl.pkl
  hde_ene.pene
  hde_en.pen
  hde_mat.pmat
  hde_dn.pdn
  hde_c.pc
  hde_cme.pcme
  hde_cm.pcm
  hde_cene.pcene
  hde_cen.pcen

* hicksian demand equations - determine demand for underlying commodities
  hde_ld_ex.ld_ex
  hde_kd_ex.kd_ex
  hde_id_ex.id_ex
  hde_res_ex.res_ex
  hde_res.res
  hde_ld.ld
  hde_kd.kd
  hde_id_ele.id
  hde_id_m.id
  hde_id_en.id
  hde_m.m
  hde_d.d
  hde_x.x
  hde_leis.leis
  hde_cd_trn.cd
  hde_cd_ele.cd
  hde_cd_m.cd
  hde_cd_en.cd
  hde_i.i
  hde_g.g

* budget constraints
  bc_hh.lambda
  bc_gov.incadj

* first order conditions
  foc_cl.cl
  foc_lambda.kh

* unit cost equations
  uce_y_ex.uc_y_ex
  uce_y.uc_y
  uce_klem.uc_klem
  uce_kle.uc_kle
  uce_kl.uc_kl
  uce_ene.uc_ene
  uce_en.uc_en
  uce_mat.uc_mat
  uce_a.uc_a
  uce_dn.uc_dn
  uce_dx.py
  uce_k_ex.pr_ex_agg
  uce_cl.uc_cl
  uce_c.uc_c
  uce_cme.uc_cme
  uce_cene.uc_cene
  uce_cen.uc_cen
  uce_cm.uc_cm
  uce_inv.uc_inv
  uce_gov.uc_gov

* closure rules
  cr_gov.pgov
  cr_kt.kt

* auxilary equations
  aux_k_ex.k_ex
  aux_w.w
  aux_cpi.cpi

/;



* ------------------------------------------------------------------------------
* set options for the solver
* ------------------------------------------------------------------------------

$label solver_options

* use the path options file
sage.optfile = 1;

* use the modeling scaling
sage.scaleopt = 1;

* iteration limit
sage.iterlim = 1000000;

* time limit
sage.reslim = 2*60*60;



* if balanced growth path starting values are requested use them
$if %balanced_start_values%==1 $goto set_using_balanced

* if a baseline file is provided use that for the starting values
$if not %gdx_baseline_file%==0 $goto set_with_saved_values



* ------------------------------------------------------------------------------
* set the starting values for the endogenous variables
* ------------------------------------------------------------------------------

$label set_using_balanced

* parameters for starting values
parameters
  q_base(t)             quantity index on balanced growth path
  p_base(t)             price index on balanced growth path
  k_ext_base(t,r)       extant capital path
  k_new_base(t,r,s)     new capital path;

* quantities grow at the steady state growth rate
q_base(t) = (1+gamma+omega)**(ord(t)-1);

* prices decline with the interest rate
p_base(t) = 1/(1+rbar)**(ord(t)-1);

* extant capital path
k_ext_base(t,r)$sum(s, extant_share(r,s)) = (1-delta_ex)**(ord(t)-1);

* new capital path
k_new_base(t,r,s)$(1-extant_share(r,s)) = ((1+gamma+omega)**(ord(t)-1)
                                           -(1-delta_ex)**(ord(t)-1)*extant_share(r,s)
                                          )/(1-extant_share(r,s));

* prices
py.l(t,r,s)          = p_base(t);
pres.l(t,r,s)        = p_base(t);
pklem.l(t,r,s)       = p_base(t);
pkle.l(t,r,s)        = p_base(t);
pkl.l(t,r,s)         = p_base(t);
pene.l(t,r,s)        = p_base(t);
pen.l(t,r,s)         = p_base(t);
pmat.l(t,r,s)        = p_base(t);
pa.l(t,r,s)          = p_base(t);
pdn.l(t,r,s)         = p_base(t);
pd.l(t,r,s)          = p_base(t);
pn.l(t,s)            = p_base(t);
pcl.l(t,r,h)         = p_base(t);
pc.l(t,r,h)          = p_base(t);
pcme.l(t,r,h)        = p_base(t);
pcene.l(t,r,h)       = p_base(t);
pcen.l(t,r,h)        = p_base(t);
pcm.l(t,r,h)         = p_base(t);
pgov.l(t,r)          = p_base(t);
pfx.l(t)             = p_base(t);
pl.l(t,r)            = p_base(t);
pk.l(t,r)            = p_base(t)/(sum(tt$tfirst(tt), p_base(tt+1))+1$(card(t) eq 1));
pr.l(t,r)            = p_base(t);
pr_ex.l(t,r,s)       = p_base(t);
pr_ex_agg.l(t,r)     = p_base(t);
pkt.l(r)             = sum(t$tlast(t), pk.l(t,r))*sum(tt$tfirst(tt), p_base(tt+1))
                       +1$(card(t) eq 1);

* shadow prices
lambda.l(t,r,h)      = (cl0(r,h)/sum(tfirst, n(tfirst,r,h)))**(-eta(r,h))/beta(t,r,h);

* activities
y_ex.l(t,r,s)        = k_ext_base(t,r);
y.l(t,r,s)           = k_new_base(t,r,s);
a.l(t,r,s)           = q_base(t);
inv.l(t,r)           = q_base(t);
gov.l(t,r)           = q_base(t);
l.l(t,r,h)           = q_base(t);

* demands
ld_ex.l(t,r,s)       = k_ext_base(t,r);
kd_ex.l(t,r,s)       = k_ext_base(t,r);
id_ex.l(t,r,ss,s)    = k_ext_base(t,r);
res_ex.l(t,r,s)      = k_ext_base(t,r);
klem.l(t,r,s)        = k_new_base(t,r,s);
res.l(t,r,s)         = k_new_base(t,r,s);
kl.l(t,r,s)          = k_new_base(t,r,s);
ld.l(t,r,s)          = k_new_base(t,r,s);
kd.l(t,r,s)          = k_new_base(t,r,s);
mat.l(t,r,s)         = k_new_base(t,r,s);
kle.l(t,r,s)         = k_new_base(t,r,s);
ene.l(t,r,s)         = k_new_base(t,r,s);
en.l(t,r,s)          = k_new_base(t,r,s);
id.l(t,r,ss,s)       = k_new_base(t,r,s);
d.l(t,r,s)           = q_base(t);
dn.l(t,r,s)          = q_base(t);
m.l(t,r,s,trd)       = q_base(t);
x.l(t,r,s,trd)       = q_base(t);
cl.l(t,r,h)          = q_base(t);
leis.l(t,r,h)        = q_base(t);
c.l(t,r,h)           = q_base(t);
cme.l(t,r,h)         = q_base(t);
cene.l(t,r,h)        = q_base(t);
cen.l(t,r,h)         = q_base(t);
cm.l(t,r,h)          = q_base(t);
cd.l(t,r,s,h)        = q_base(t);
i.l(t,r,s)           = q_base(t);
g.l(t,r,s)           = q_base(t);

* stocks
k_ex.l(t,r)          = k_ext_base(t,r);
k.l(t,r)             = sum(s, kd.l(t,r,s)*kd0(r,s))/r0(r);
kh.l(t,r,h)          = sum(s, kd.l(t,r,s)*kd0(r,s))/r0(r);
kt.l(r)              = sum(s$(1-extant_share(r,s)), ((1+gamma+omega)**card(t)
                                                     -(1-delta_ex)**card(t)*extant_share(r,s)
                                                    )/(1-extant_share(r,s))*kd0(r,s)
                          )/r0(r);

* unit costs
uc_y_ex.l(t,r,s)     = (1-ty0(r,s))*p_base(t);
uc_y.l(t,r,s)        = (1-ty0(r,s))*p_base(t);
uc_klem.l(t,r,s)     = p_base(t);
uc_kl.l(t,r,s)       = p_base(t);
uc_mat.l(t,r,s)      = p_base(t);
uc_kle.l(t,r,s)      = p_base(t);
uc_ene.l(t,r,s)      = p_base(t);
uc_en.l(t,r,s)       = p_base(t);
uc_a.l(t,r,s)        = p_base(t);
uc_dn.l(t,r,s)       = p_base(t);
uc_cl.l(t,r,h)       = p_base(t);
uc_c.l(t,r,h)        = p_base(t);
uc_cme.l(t,r,h)      = p_base(t);
uc_cene.l(t,r,h)     = p_base(t);
uc_cen.l(t,r,h)      = p_base(t);
uc_cm.l(t,r,h)       = p_base(t);
uc_inv.l(t,r)        = p_base(t);
uc_gov.l(t,r)        = p_base(t);

* auxillary
incadj.l(t)          = 0;
w.l(r,h)             = sum(t,  n(t,r,h)*beta(t,r,h)*(cl.l(t,r,h)*cl0(r,h)
                                                      /n(t,r,h)
                                                    )**(1-eta(r,h))/(1-eta(r,h)));
cpi.l(t)             = p_base(t);

$goto scale_equations



* ------------------------------------------------------------------------------
* set the starting values for the endogenous variables using baseline results
* ------------------------------------------------------------------------------

$label set_with_saved_values

* open the gdx file connection
$gdxin %gdx_baseline_file%

* prices
$loadm py
$loadm pres
$loadm pklem
$loadm pkle
$loadm pkl
$loadm pene
$loadm pen
$loadm pmat
$loadm pa
$loadm pdn
$loadm pd
$loadm pn
$loadm pcl
$loadm pc
$loadm pcme
$loadm pcene
$loadm pcen
$loadm pcm
$loadm pgov
$loadm pfx
$loadm pl
$loadm pk
$loadm pr
$loadm pr_ex
$loadm pr_ex_agg
$loadm pkt

* shadow prices
$loadm lambda

* activities
$loadm y_ex
$loadm y
$loadm a
$loadm inv
$loadm gov
$loadm l

* demands
$loadm ld_ex
$loadm kd_ex
$loadm id_ex
$loadm res_ex
$loadm klem
$loadm res
$loadm kl
$loadm ld
$loadm kd
$loadm mat
$loadm kle
$loadm ene
$loadm en
$loadm id
$loadm d
$loadm dn
$loadm m
$loadm x
$loadm cl
$loadm leis
$loadm c
$loadm cme
$loadm cene
$loadm cen
$loadm cm
$loadm cd
$loadm i
$loadm g

* stocks
$loadm k_ex
$loadm k
$loadm kh
$loadm kt

* unit costs
$loadm uc_y_ex
$loadm uc_y
$loadm uc_klem
$loadm uc_kl
$loadm uc_mat
$loadm uc_kle
$loadm uc_ene
$loadm uc_en
$loadm uc_a
$loadm uc_dn
$loadm uc_cl
$loadm uc_c
$loadm uc_cme
$loadm uc_cene
$loadm uc_cen
$loadm uc_cm
$loadm uc_inv
$loadm uc_gov

* auxillary
$loadm incadj
$loadm w
$loadm cpi

* make sure none of the prices are zero which can cause a division by zero error
* prices can end up at zero if the model zeros out production in a region,
* sector, capital vintage combination
py.l(t,r,s)$(not py.l(t,r,s))           = 1e-6;
pres.l(t,r,s)$(not pres.l(t,r,s))       = 1e-6;
pklem.l(t,r,s)$(not pklem.l(t,r,s))     = 1e-6;
pkl.l(t,r,s)$(not pkl.l(t,r,s))         = 1e-6;
pmat.l(t,r,s)$(not pmat.l(t,r,s))       = 1e-6;
pkle.l(t,r,s)$(not pkle.l(t,r,s))       = 1e-6;
pene.l(t,r,s)$(not pene.l(t,r,s))       = 1e-6;
pen.l(t,r,s)$(not pen.l(t,r,s))         = 1e-6;
pr_ex.l(t,r,s)$(not pr_ex.l(t,r,s))     = 1e-6;
pr_ex_agg.l(t,r)$(not pr_ex_agg.l(t,r)) = 1e-6;



* ------------------------------------------------------------------------------
* scale the model
* ------------------------------------------------------------------------------

$label scale_equations

bc_hh.scale(t,r,h)           = 1e4;
mc_a.scale(t,r,s)            = 1e3;
zp_k.scale(t,r)              = 1e3;
bc_gov.scale(t)              = 1e3;
mc_pk.scale(t,r)             = 1e3;
mc_pkt.scale(r)              = 1e3;
mc_l.scale(t,r)              = 1e3;
mc_pr_ex.scale(t,r,s)        = 1e1;
mc_pr.scale(t,r)             = 1e3;
mc_time.scale(t,r,h)         = 1e3;
mc_fx.scale(t)               = 1e3;
cr_gov.scale(t,r)            = 1e2;
hde_d.scale(t,r,s)           = 1e3;
hde_dn.scale(t,r,s)          = 1e2;
hde_m.scale(t,r,s,trd)       = 1e3;
hde_x.scale(t,r,s,trd)       = 1e2;
hde_kd.scale(t,r,s)          = 1e1;
hde_en.scale(t,r,s)          = 1e1;
hde_kl.scale(t,r,s)          = 1e1;
hde_kle.scale(t,r,s)         = 1e1;
hde_mat.scale(t,r,s)         = 1e1;
hde_ld.scale(t,r,s)          = 1e1;
mc_dtrd.scale(t,s)           = 1e2;
mc_d.scale(t,r,s)            = 1e2;
hde_leis.scale(t,r,h)        = 1e2;
hde_c.scale(t,r,h)           = 1e1;
hde_cme.scale(t,r,h)         = 1e1;
hde_cene.scale(t,r,h)        = 1e1;
hde_cen.scale(t,r,h)         = 1e1;
hde_cm.scale(t,r,h)          = 1e1;
hde_cd_trn.scale(t,r,stn,h)  = 1e1;
hde_cd_en.scale(t,r,scen,h)  = 1e1;
hde_cd_ele.scale(t,r,sel,h)  = 1e1;
hde_cd_m.scale(t,r,scm,h)    = 1e1;
hde_res.scale(t,r,s)         = 1e1;
hde_ene.scale(t,r,s)         = 1e2;
foc_cl.scale(t,r,h)          = 1e2;



* ------------------------------------------------------------------------------
* fixed values in most applications and other starting value issues
* ------------------------------------------------------------------------------

* set the numeraire
pfx.fx(t)$tfirst(t) = 1;

* initial household capital stock is fixed
kh.fx(t,r,h)$tfirst(t) = 1;

* small open economy model
pm.fx(t,s) = 1;
px.fx(t,s) = 1;

* allow for the option to perturb the starting values (useful for testing)
y.l(t,r,s) = y.l(t,r,s)+%perturb_start%;



* ------------------------------------------------------------------------------
* if this is a static run fix additional variables
* ------------------------------------------------------------------------------

if (card(t) eq 1,

* fix capital since its not used
  k.fx(t,r) = 1;

* fix terminal capital so that real investment is fixed
  kt.fx(r) = inv0(r)/k0(r);

);



* ------------------------------------------------------------------------------
* write solver options
* ------------------------------------------------------------------------------

* open the path options file for writing
file path_opt /path.opt/;
put path_opt;

* write the path options
put /"convergence_tolerance 1e-10;"/;
put /"output_final_scaling_statistics yes;"/;
put /"nms_initial_reference_factor 2;"/;

* turning off the crash phase is usually the fastest way to solve the model
put /"crash_method none;"/;

* close the path options file
putclose;



* ------------------------------------------------------------------------------
* set any policy shocks
* ------------------------------------------------------------------------------

$if not %policy_file%==0 $include %policy_file%



* ------------------------------------------------------------------------------
* solve the model
* ------------------------------------------------------------------------------

* solve the model
solve sage using mcp;



* ------------------------------------------------------------------------------
* save the results
* ------------------------------------------------------------------------------

if (%gdx_save% eq 1, execute_unload "%gdx_results_file%");



* ------------------------------------------------------------------------------
* write output files with the results
* ------------------------------------------------------------------------------

* additional output variables
parameter
    gdp(t)              real gross dometic product
    incadjtot(t,r,h)    full lump sum transfer
    tax_revenue(t,r)    total tax revenue;

* real gdp
gdp(t) = sum(tt$tfirst(tt),
               sum((r,s,h), pa.l(tt,r,s)*cd.l(t,r,s,h)*cd_base(t,r,s,h))
               +sum((r,s), pa.l(tt,r,s)*i.l(t,r,s)*i0(r,s))
               +sum((r,s), pa.l(tt,r,s)*g.l(t,r,s)*g_base(t,r,s))
               +sum((r,s), pfx.l(tt)*px.l(tt,s)*x.l(t,r,s,"ftrd")*x0(r,s,"ftrd"))
               -sum((r,s), pfx.l(tt)*pm.l(tt,s)*m.l(t,r,s,"ftrd")*m0(r,s,"ftrd")));

* income adjustment - lump sum tax/transfer
incadjtot(t,r,h) =  cpi.l(t)*incadj0(r,h)*(1+gamma+omega)**(ord(t)-1)
                   +cpi.l(t)*incadj.l(t)*incadj_share(r,h);

* total tax tax revenue
tax_revenue(t,r) =  sum(s, sum(h, tc(t,r)*pa.l(t,r,s)*cd.l(t,r,s,h)*cd_base(t,r,s,h))
                          +ty_ex(t,r,s)*py.l(t,r,s)*y_ex.l(t,r,s)*y_ex0(r,s)
                          +ty(t,r,s)*py.l(t,r,s)*y.l(t,r,s)*y0(r,s)
                          +tl(t,r)*pl.l(t,r)*ld_ex.l(t,r,s)*ld_ex0(r,s)
                          +tl(t,r)*pl.l(t,r)*ld.l(t,r,s)*ld_base(t,r,s)
                          +tk(t,r)*pr_ex.l(t,r,s)*kd_ex.l(t,r,s)*kd_ex0(r,s)
                          +tk(t,r)*pr.l(t,r)*kd.l(t,r,s)*kd_base(t,r,s)
                          +tk(t,r)*pres.l(t,r,s)*( res_ex.l(t,r,s)*res_ex0(r,s)
                                                  +res.l(t,r,s)*res0(r,s)));

* output csv results file
file results /%output_file%/;
results.pc = 5;
results.nd = 6;
put results;

put "parameter", "time", "region", "sector", "household", "sector.from", "value"/;

* relative prices
loop((t,r,s),   put   "py",       t.tl,  r.tl,  s.tl,  "",    "",    py.l(t,r,s)       / );
loop((t,r,s),   put   "pa",       t.tl,  r.tl,  s.tl,  "",    "",    pa.l(t,r,s)       / );
loop((t,r,h),   put   "pc",       t.tl,  r.tl,  "",    h.tl,  "",    pc.l(t,r,h)       / );
loop((t,r,h),   put   "pcl",      t.tl,  r.tl,  "",    h.tl,  "",    pcl.l(t,r,h)      / );
loop((t,r),     put   "pl",       t.tl,  r.tl,  "",    "",    "",    pl.l(t,r)         / );
loop((t,r,s),   put   "pr_ex",    t.tl,  r.tl,  s.tl,  "",    "",    pr_ex.l(t,r,s)    / );
loop((t,r,s),   put   "pres",     t.tl,  r.tl,  s.tl,  "",    "",    pres.l(t,r,s)     / );
loop((t,r),     put   "pr",       t.tl,  r.tl,  "",    "",    "",    pr.l(t,r)         / );
loop((t,r),     put   "pk",       t.tl,  r.tl,  "",    "",    "",    pk.l(t,r)         / );
loop((t,r),     put   "pgov",     t.tl,  r.tl,  "",    "",    "",    pgov.l(t,r)       / );
loop((t,s),     put   "pn",       t.tl,  "",    s.tl,  "",    "",    pn.l(t,s)         / );
loop((t,r,s),   put   "pd",       t.tl,  r.tl,  s.tl,  "",    "",    pd.l(t,r,s)       / );
loop((t),       put   "pfx",      t.tl,  "",    "",    "",    "",    pfx.l(t)          / );
loop((r),       put   "pkt",      "",    r.tl,  "",    "",    "",    pkt.l(r)          / );
loop((t),       put   "cpi",      t.tl,  "",    "",    "",    "",    cpi.l(t)          / );

* quantities
loop((t,r,s),   put   "y_ex",     t.tl,  r.tl,  s.tl,  "",     "",    (y_ex.l(t,r,s)*y_ex0(r,s))         / );
loop((t,r,s),   put   "y",        t.tl,  r.tl,  s.tl,  "",     "",    (y.l(t,r,s)*y0(r,s))               / );
loop((t,r,s),   put   "a",        t.tl,  r.tl,  s.tl,  "",     "",    (a.l(t,r,s)*a0(r,s))               / );
loop((t,r,s),   put   "d",        t.tl,  r.tl,  s.tl,  "",     "",    (d.l(t,r,s)*d0(r,s))               / );
loop((t,r,s),   put   "i",        t.tl,  r.tl,  s.tl,  "",     "",    (i.l(t,r,s)*i0(r,s))               / );
loop((t,r,s,h), put   "cd",       t.tl,  r.tl,  s.tl,  h.tl,   "",    (cd.l(t,r,s,h)*cd_base(t,r,s,h))   / );
loop((t,r,h),   put   "c",        t.tl,  r.tl,  "",    h.tl,   "",    (c.l(t,r,h)*c0(r,h))               / );
loop((t,r,h),   put   "cl",       t.tl,  r.tl,  "",    h.tl,   "",    (cl.l(t,r,h)*cl0(r,h))             / );
loop((t,r,s),   put   "g",        t.tl,  r.tl,  s.tl,  "",     "",    (g.l(t,r,s)*g_base(t,r,s))         / );
loop((t,r,s),   put   "x_f",      t.tl,  r.tl,  s.tl,  "",     "",    (x.l(t,r,s,"ftrd")*x0(r,s,"ftrd")) / );
loop((t,r,s),   put   "m_f",      t.tl,  r.tl,  s.tl,  "",     "",    (m.l(t,r,s,"ftrd")*m0(r,s,"ftrd")) / );
loop((t,r,s),   put   "x_d",      t.tl,  r.tl,  s.tl,  "",     "",    (x.l(t,r,s,"dtrd")*x0(r,s,"dtrd")) / );
loop((t,r,s),   put   "m_d",      t.tl,  r.tl,  s.tl,  "",     "",    (m.l(t,r,s,"dtrd")*m0(r,s,"dtrd")) / );
loop((t,r,s),   put   "ld_ex",    t.tl,  r.tl,  s.tl,  "",     "",    (ld_ex.l(t,r,s)*ld_ex0(r,s))       / );
loop((t,r,s),   put   "ld",       t.tl,  r.tl,  s.tl,  "",     "",    (ld.l(t,r,s)*ld_base(t,r,s))       / );
loop((t,r,h),   put   "leisure",  t.tl,  r.tl,  "",    h.tl,   "",    (leis.l(t,r,h)*leis0(r,h))         / );
loop((t,r,h),   put   "l",        t.tl,  r.tl,  "",    h.tl,   "",    (l.l(t,r,h)*l0(r,h))               / );
loop((t,r,s),   put   "kd_ex",    t.tl,  r.tl,  s.tl,  "",     "",    (kd_ex.l(t,r,s)*kd_ex0(r,s))       / );
loop((t,r,s),   put   "kd",       t.tl,  r.tl,  s.tl,  "",     "",    (kd.l(t,r,s)*kd_base(t,r,s))       / );
loop((t,r,h),   put   "r_ex",     t.tl,  r.tl,  "",    h.tl,   "",    (k_ex.l(t,r)*re_ex0(r,h))          / );
loop((t,r,s,h), put   "res_ex",   t.tl,  r.tl,  s.tl,  h.tl,   "",    (res_ex.l(t,r,s)*rese_ex0(r,s,h))  / );
loop((t,r,s,h), put   "res",      t.tl,  r.tl,  s.tl,  h.tl,   "",    (res.l(t,r,s)*rese0(r,s,h))        / );

loop((t,r,s,ss),put   "id",       t.tl,  r.tl,  s.tl,  "",   ss.tl,   (id.l(t,r,ss,s)*id_base(t,r,ss,s)) / );
loop((t,r,s,ss),put   "id_ex",    t.tl,  r.tl,  s.tl,  "",   ss.tl,   (id_ex.l(t,r,ss,s)*id_ex0(r,ss,s)) / );


loop((t,r),     put   "k",        t.tl,  r.tl,  "",    "",     "",    (k.l(t,r)*k0(r))                   / );
loop((r),       put   "kt",       "",    r.tl,  "",    "",     "",    (kt.l(r)*k0(r))                    / );
loop((t,r,h),   put   "kh",       t.tl,  r.tl,  "",    h.tl,   "",    (kh.l(t,r,h)*ke0(r,h))             / );

loop((r,h),     put   "w",        "",    r.tl,  "",    h.tl,   "",    w.l(r,h)                           / );


* additional output
loop((t),       put   "gdp",      t.tl,  "",    "",    "",     "",    gdp(t)              / );
loop((t,r,h),   put   "incadj",   t.tl,  r.tl,  "",    h.tl,   "",    incadjtot(t,r,h)    / );
loop((t,r,h),   put   "bopdef",   t.tl,  r.tl,  "",    h.tl,   "",    (bopdef(r,h)*n(t,r,h)/sum(tt$tfirst(tt), n(tt,r,h)))    / );
loop((t,r,h),   put   "te",       t.tl,  r.tl,  "",    h.tl,   "",    te(t,r,h)           / );
loop((t,r,h),   put   "n",        t.tl,  r.tl,  "",    h.tl,   "",    n(t,r,h)           / );

* taxes
loop((t,r),     put   "tk",          t.tl,  r.tl,  "",    "",    "",     tk(t,r)          / );
loop((t,r),     put   "tl",          t.tl,  r.tl,  "",    "",    "",     tl(t,r)          / );
loop((t,r),     put   "tc",          t.tl,  r.tl,  "",    "",    "",     tc(t,r)          / );
loop((t,r,s),   put   "ty",          t.tl,  r.tl,  s.tl,  "",    "",     ty(t,r,s)        / );
loop((t,r),     put   "tax_revenue", t.tl,  r.tl,  "",    "",    "",     tax_revenue(t,r) / );

putclose;



* ------------------------------------------------------------------------------
* display additional statistics for the model run
* ------------------------------------------------------------------------------

parameters
  inv0_end(r)           endogenous benchmark investment
  c0_end(r,h)           endogenous benchmark consumption;

inv0_end(r) = sum((t,s)$tfirst(t), pa.l(t,r,s)*i.l(t,r,s)*i0(r,s));
display inv0, inv0_end;

c0_end(r,h) = sum((t,s)$tfirst(t), (1+tc(t,r))*pa.l(t,r,s)*cd.l(t,r,s,h)
                                                          *cd_base(t,r,s,h));
display c0, c0_end;

* calculate the implicit annual interest rate in each region
parameters
  r_annual(r,h)         calibrated annual interest rate;
r_annual(r,h) = (1+r0(r)/k0(r)-delta)**(1/interval)-1;
display r_annual, k0;

* display gdp growth
parameters
  gdp_growth(t)         endogenous gdp growth;
gdp_growth(t)$(ord(t) gt 1) = (gdp(t)/gdp(t-1))**(1/interval)-1;
display gdp_growth;

* present value costs of transforming extant capital
parameters
  k_ex_cet_costs(t,r)   present value costs of transforming extant capital;
k_ex_cet_costs(t,r) = pr_ex_agg.l(t,r)*(k_ex.l(t,r)*k_ex0(r)
                                        -sum(s, kd_ex.l(t,r,s)*kd_ex0(r,s)));
display k_ex_cet_costs;



* ------------------------------------------------------------------------------
* run solution checks
* ------------------------------------------------------------------------------

parameters
  zero_tol         tolerance for assuming a zero difference        / 1e-4 /;



* ------------------------------------------------------------------------------
* gdp calculation
* ------------------------------------------------------------------------------

* Calculation of nominal GDP from the expenditure and value added sides
parameters
  gdp_e(t)           nominal gdp - calculated based on expenditures
  gdp_va(t)          nominal gdp - calculated based on value added;

* nominal gdp - calculated based on expenditures
gdp_e(t)  = sum(r, sum((h,s), pa.l(t,r,s)*(cd.l(t,r,s,h)*cd_base(t,r,s,h)))
            +sum(s, pa.l(t,r,s)*(g.l(t,r,s)*g_base(t,r,s)+i.l(t,r,s)*i0(r,s)))
            +sum((s,trd), pfx.l(t)*( px.l(t,s)*x.l(t,r,s,trd)*x0(r,s,trd)
                                    -pm.l(t,s)*m.l(t,r,s,trd)*m0(r,s,trd))));

* nominal gdp - calculated based on value added
gdp_va(t) = sum(r, sum(s, pl.l(t,r)*ld_ex.l(t,r,s)*ld_ex0(r,s)
                         +pl.l(t,r)*ld.l(t,r,s)*ld_base(t,r,s)
                         +pr_ex.l(t,r,s)*kd_ex.l(t,r,s)*kd_ex0(r,s)
                         +pr.l(t,r)*kd.l(t,r,s)*kd_base(t,r,s)
                         +pres.l(t,r,s)*(res_ex.l(t,r,s)*res_ex0(r,s)
                                         +res.l(t,r,s)*res0(r,s))
                      )
                   +sum(s, tl(t,r)*pl.l(t,r)*ld_ex.l(t,r,s)*ld_ex0(r,s)
                          +tl(t,r)*pl.l(t,r)*ld.l(t,r,s)*ld_base(t,r,s)
                          +tk(t,r)*pr_ex.l(t,r,s)*kd_ex.l(t,r,s)*kd_ex0(r,s)
                          +tk(t,r)*pr.l(t,r)*kd.l(t,r,s)*kd_base(t,r,s)
                          +tk(t,r)*pres.l(t,r,s)*(res_ex.l(t,r,s)*res_ex0(r,s)
                                                  +res.l(t,r,s)*res0(r,s))
                          +ty_ex(t,r,s)*py.l(t,r,s)*y_ex.l(t,r,s)*y_ex0(r,s)
                          +ty(t,r,s)*py.l(t,r,s)*y.l(t,r,s)*y0(r,s)));

 display gdp_e, gdp_va;

* savings rate
parameters
  savings_rate(t);
savings_rate(t) = sum((r,s), pa.l(t,r,s)*i.l(t,r,s)*i0(r,s))/gdp_e(t);
display savings_rate;



* ------------------------------------------------------------------------------
* sam - build the social accounting matrix for each region
* ------------------------------------------------------------------------------

parameter
  sam_row_comm(t,r,*)   sam - row - commodity account
  sam_col_comm(t,r,*)   sam - col - commodity account
  sam_row_act(t,r,*)    sam - row - activity account
  sam_col_act(t,r,*)    sam - col - activity account
  sam_row_govt(t,*)     sam - row - government account
  sam_col_govt(t,*)     sam - col - government account
  sam_row_hh(t,r,*)     sam - row - household account
  sam_col_hh(t,r,*)     sam - col - household account
  sam_row_row(t,r,*)    sam - row - rest of world account
  sam_col_row(t,r,*)    sam - col - rest of world account
  sam_table(t,*,*)      sam - summary table;

* sam - row - commodity account
sam_row_comm(t,r,"pa*id")        = sum((ss,s), pa.l(t,r,s)*id_ex.l(t,r,s,ss)*id_ex0(r,s,ss)
                                              +pa.l(t,r,s)*id.l(t,r,s,ss)*id_base(t,r,s,ss));
sam_row_comm(t,r,"pa*cd")        = sum((s,h), pa.l(t,r,s)*cd.l(t,r,s,h)*cd_base(t,r,s,h));
sam_row_comm(t,r,"pa*g")         = sum(s, pa.l(t,r,s)*g.l(t,r,s)*g_base(t,r,s));
sam_row_comm(t,r,"pa*i")         = sum(s, pa.l(t,r,s)*i.l(t,r,s)*i0(r,s));
sam_row_comm(t,r,"comm row tot") =  sam_row_comm(t,r,"pa*id")
                                   +sam_row_comm(t,r,"pa*cd")
                                   +sam_row_comm(t,r,"pa*g")
                                   +sam_row_comm(t,r,"pa*i");

* sam - col - commodity account
sam_col_comm(t,r,"pd*d")         = sum(s, py.l(t,r,s)*(y.l(t,r,s)*y0(r,s)
                                                       +y_ex.l(t,r,s)*y_ex0(r,s)))
                                   -sum((s,ftrd), pfx.l(t)*px.l(t,s)*x.l(t,r,s,ftrd)*x0(r,s,ftrd))
                                   -sum((s,dtrd), pn.l(t,s)*x.l(t,r,s,dtrd)*x0(r,s,dtrd));
sam_col_comm(t,r,"pfx*m_f")      = sum(s, pfx.l(t)*pm.l(t,s)*m.l(t,r,s,"ftrd")*m0(r,s,"ftrd"));
sam_col_comm(t,r,"pn*m_n")       = sum(s, pn.l(t,s)*m.l(t,r,s,"dtrd")*m0(r,s,"dtrd"));
sam_col_comm(t,r,"comm col tot") =  sam_col_comm(t,r,"pd*d")
                                   +sam_col_comm(t,r,"pfx*m_f")
                                   +sam_col_comm(t,r,"pn*m_n");
sam_col_comm(t,r,"comm dif")     =  sam_row_comm(t,r,"comm row tot")
                                   -sam_col_comm(t,r,"comm col tot");

* sam - row - activity account
sam_row_act(t,r,"pd*d")         = sum(s, py.l(t,r,s)*(y.l(t,r,s)*y0(r,s)
                                                      +y_ex.l(t,r,s)*y_ex0(r,s)))
                                  -sum((s,ftrd), pfx.l(t)*px.l(t,s)*x.l(t,r,s,ftrd)*x0(r,s,ftrd))
                                  -sum((s,dtrd), pn.l(t,s)*x.l(t,r,s,dtrd)*x0(r,s,dtrd));
sam_row_act(t,r,"pfx*x_f")      =  sum(s, pfx.l(t)*px.l(t,s)*x.l(t,r,s,"ftrd")*x0(r,s,"ftrd"));
sam_row_act(t,r,"pn*x_n")       = sum(s, pn.l(t,s)*x.l(t,r,s,"dtrd")*x0(r,s,"dtrd"));
sam_row_act(t,r,"act row tot")  =  sam_row_act(t,r,"pd*d")
                                  +sam_row_act(t,r,"pfx*x_f")
                                  +sam_row_act(t,r,"pn*x_n");

* sam - col - activity account
sam_col_act(t,r,"pa*id")       = sum((ss,s), pa.l(t,r,s)*id_ex.l(t,r,s,ss)*id_ex0(r,s,ss)
                                            +pa.l(t,r,s)*id.l(t,r,s,ss)*id_base(t,r,s,ss));
sam_col_act(t,r,"pl*ld")       = sum(s, pl.l(t,r)*ld_ex.l(t,r,s)*ld_ex0(r,s)
                                       +pl.l(t,r)*ld.l(t,r,s)*ld_base(t,r,s));
sam_col_act(t,r,"pr*kd")       = sum(s, pr_ex.l(t,r,s)*kd_ex.l(t,r,s)*kd_ex0(r,s)
                                       +pr.l(t,r)*kd.l(t,r,s)*kd_base(t,r,s)
                                       +pres.l(t,r,s)*(res_ex.l(t,r,s)*res_ex0(r,s)
                                                       +res.l(t,r,s)*res0(r,s)));
sam_col_act(t,r,"tax")         =  sum(s, ty_ex(t,r,s)*py.l(t,r,s)*y_ex.l(t,r,s)*y_ex0(r,s)
                                        +ty(t,r,s)*py.l(t,r,s)*y.l(t,r,s)*y0(r,s))
                                 +sum(s, tl(t,r)*pl.l(t,r)*ld_ex.l(t,r,s)*ld_ex0(r,s)
                                        +tl(t,r)*pl.l(t,r)*ld.l(t,r,s)*ld_base(t,r,s))
                                 +sum(s, tk(t,r)*pr_ex.l(t,r,s)*kd_ex.l(t,r,s)*kd_ex0(r,s)
                                        +tk(t,r)*pr.l(t,r)*kd.l(t,r,s)*kd_base(t,r,s)
                                        +tk(t,r)*pres.l(t,r,s)*(res_ex.l(t,r,s)*res_ex0(r,s)
                                                                +res.l(t,r,s)*res0(r,s)));
sam_col_act(t,r,"act col tot") =  sam_col_act(t,r,"pa*id")
                                 +sam_col_act(t,r,"pl*ld")
                                 +sam_col_act(t,r,"pr*kd")
                                 +sam_col_act(t,r,"tax");
sam_col_act(t,r,"act dif")     =  sam_row_act(t,r,"act row tot")
                                 -sam_col_act(t,r,"act col tot");

* sam - row - government account
sam_row_govt(t,"c tax")        = sum((r,h,s), tc(t,r)*pa.l(t,r,s)*cd.l(t,r,s,h)*cd_base(t,r,s,h));
sam_row_govt(t,"y tax")        = sum((r,s), ty_ex(t,r,s)*py.l(t,r,s)*y_ex.l(t,r,s)*y_ex0(r,s)
                                           +ty(t,r,s)*py.l(t,r,s)*y.l(t,r,s)*y0(r,s));
sam_row_govt(t,"l tax")        = sum((r,s), tl(t,r)*pl.l(t,r)*ld_ex.l(t,r,s)*ld_ex0(r,s)
                                           +tl(t,r)*pl.l(t,r)*ld.l(t,r,s)*ld_base(t,r,s));
sam_row_govt(t,"k tax")        = sum((r,s), tk(t,r)*pr_ex.l(t,r,s)*kd_ex.l(t,r,s)*kd_ex0(r,s)
                                           +tk(t,r)*pr.l(t,r)*kd.l(t,r,s)*kd_base(t,r,s)
                                           +tk(t,r)*pres.l(t,r,s)*(res_ex.l(t,r,s)*res_ex0(r,s)
                                                               +res.l(t,r,s)*res0(r,s)));
sam_row_govt(t,"govt row tot") =  sam_row_govt(t,"c tax")
                                 +sam_row_govt(t,"y tax")
                                 +sam_row_govt(t,"l tax")
                                 +sam_row_govt(t,"k tax");

* sam - col - government account
sam_col_govt(t,"gov purchases") = sum((r,s), pa.l(t,r,s)*g_base(t,r,s));
sam_col_govt(t,"gov purchases") = sum(r, pgov.l(t,r)*gov.l(t,r)*gov0(r));
sam_col_govt(t,"gov exo trans") = sum((r,h), cpi.l(t)*incadj0(r,h)*(1+gamma+omega)**(ord(t)-1));
sam_col_govt(t,"gov end trans") = cpi.l(t)*incadj.l(t);
sam_col_govt(t,"govt col tot")  =  sam_col_govt(t,"gov purchases")
                                  +sam_col_govt(t,"gov end trans")
                                  +sam_col_govt(t,"gov exo trans");
sam_col_govt(t,"govt dif")      =  sam_row_govt(t,"govt row tot")
                                  -sam_col_govt(t,"govt col tot");

* sam - row - household account
sam_row_hh(t,r,"l inc")         = sum(h, pl.l(t,r)*l.l(t,r,h)*l0(r,h));
sam_row_hh(t,r,"k inc")         = pr_ex_agg.l(t,r)*k_ex.l(t,r)*k_ex0(r)
                                  +sum(s, pr.l(t,r)*kd.l(t,r,s)*kd_base(t,r,s)
                                         +pres.l(t,r,s)*(res_ex.l(t,r,s)*res_ex0(r,s)
                                                         +res.l(t,r,s)*res0(r,s)));
sam_row_hh(t,r,"gov exo trans") = sum(h, cpi.l(t)*incadj0(r,h)*(1+gamma+omega)**(ord(t)-1));
sam_row_hh(t,r,"gov end trans") = sum(h, cpi.l(t)*incadj.l(t)*incadj_share(r,h));
sam_row_hh(t,r,"bopdef")        = sum(h, pfx.l(t)*bopdef(r,h)*(1+gamma+omega)**(ord(t)-1));
sam_row_hh(t,r,"hh tot")        =  sam_row_hh(t,r,"l inc")
                                  +sam_row_hh(t,r,"k inc")
                                  +sam_row_hh(t,r,"gov exo trans")
                                  +sam_row_hh(t,r,"gov end trans")
                                  +sam_row_hh(t,r,"bopdef");

* sam - col - household account
sam_col_hh(t,r,"hh cons") = sum(h, pc.l(t,r,h)*c.l(t,r,h)*c0(r,h));
sam_col_hh(t,r,"save")    = sum(s, pa.l(t,r,s)*i.l(t,r,s)*i0(r,s));
sam_col_hh(t,r,"hh tot")  =  sam_col_hh(t,r,"hh cons")
                            +sam_col_hh(t,r,"save");
sam_col_hh(t,r,"hh dif")  =  sam_row_hh(t,r,"hh tot")
                            -sam_col_hh(t,r,"hh tot");

* sam - row - rest of world account
sam_row_row(t,r,"imports_dtrd") = sum(s, pn.l(t,s)*m.l(t,r,s,"dtrd")*m0(r,s,"dtrd"));
sam_row_row(t,r,"imports_ftrd") = sum(s, pfx.l(t)*pm.l(t,s)*m.l(t,r,s,"ftrd")*m0(r,s,"ftrd"));
sam_row_row(t,r,"row tot")      =  sam_row_row(t,r,"imports_dtrd")
                                  +sam_row_row(t,r,"imports_ftrd");

* sam - col - rest of world account
sam_col_row(t,r,"exports_dtrd") = sum(s, pn.l(t,s)*x.l(t,r,s,"dtrd")*x0(r,s,"dtrd"));
sam_col_row(t,r,"exports_ftrd") = sum(s, pfx.l(t)*px.l(t,s)*x.l(t,r,s,"ftrd")*x0(r,s,"ftrd"));
sam_col_row(t,r,"bopdef")       = sum(h, pfx.l(t)*bopdef(r,h)*(1+gamma+omega)**(ord(t)-1));
sam_col_row(t,r,"row tot")      =  sam_col_row(t,r,"exports_dtrd")
                                  +sam_col_row(t,r,"exports_ftrd")
                                  +sam_col_row(t,r,"bopdef");
sam_col_row(t,r,"row dif")      =  sam_row_row(t,r,"row tot")
                                  -sam_col_row(t,r,"row tot");

* sam summary table
sam_table(t,"Commodity","Row Tot")     = sum(r, sam_row_comm(t,r,"comm row tot"));
sam_table(t,"Commodity","Col Tot")     = sum(r, sam_col_comm(t,r,"comm col tot"));
sam_table(t,"Commodity","Diff")        = sum(r, sam_col_comm(t,r,"comm dif"));
sam_table(t,"Activity","Row Tot")      = sum(r, sam_row_act(t,r,"act row tot"));
sam_table(t,"Activity","Col Tot")      = sum(r, sam_col_act(t,r,"act col tot"));
sam_table(t,"Activity","Diff")         = sum(r, sam_col_act(t,r,"act dif"));
sam_table(t,"Household","Row Tot")     = sum(r, sam_row_hh(t,r,"hh tot"));
sam_table(t,"Household","Col Tot")     = sum(r, sam_col_hh(t,r,"hh tot"));
sam_table(t,"Household","Diff")        = sum(r, sam_col_hh(t,r,"hh dif"));
sam_table(t,"Government","Row Tot")    = sam_row_govt(t,"govt row tot");
sam_table(t,"Government","Col Tot")    = sam_col_govt(t,"govt col tot");
sam_table(t,"Government","Diff")       = sam_col_govt(t,"govt dif");
sam_table(t,"Rest of World","Row Tot") = sum(r, sam_row_row(t,r,"row tot"));
sam_table(t,"Rest of World","Col Tot") = sum(r, sam_col_row(t,r,"row tot"));
sam_table(t,"Rest of World","Diff")    = sum(r, sam_col_row(t,r,"row dif"));

* output the sam
display sam_row_comm, sam_col_comm;
display sam_row_act, sam_col_act;
display sam_row_govt, sam_col_govt;
display sam_row_hh, sam_col_hh;
display sam_row_row, sam_col_row;
display sam_table;



* ------------------------------------------------------------------------------
* check for discrepency between household and market levels of capital
* ------------------------------------------------------------------------------

parameters
  k_check(t)            aggregate difference between household and market capital stock;

* national level difference between the capital stock held by households and the
* the level installed in markets - new capital is mobile across regions so the
* regional level of capital may not match that region's household capital
* holdings but at the national level these should match
k_check(t) = sum(r,pk.l(t,r)*k.l(t,r)*k0(r)
                   -sum(h, pk.l(t,r)*kh.l(t,r,h)*ke0(r,h)));
display k_check;



* ------------------------------------------------------------------------------
* output results of the solution check
* ------------------------------------------------------------------------------

set
  tnow(t)               current t in the loop;

display "----------------------------------------"
display "Solution Check Results"
display "----------------------------------------"
loop(t,
  tnow(tt) = no;
  tnow(t) = yes;
  display tnow;
  if (abs(gdp_e(t)-gdp_va(t))>zero_tol,
    display "GDP check                     FAILED";
  else
    display "GDP check                     PASSED";);
  if (abs(sam_table(t,"Commodity","Diff"))>zero_tol,
    display "Commodity account check       FAILED";
  else
    display "Commodity account check       PASSED";);
  if (abs(sam_table(t,"Activity","Diff"))>zero_tol,
    display "Activity account check        FAILED";
  else
    display "Activity account check        PASSED";);
  if (abs(sam_table(t,"Household","Diff"))>zero_tol,
    display "Household account check       FAILED";
  else
    display "Household account check       PASSED";);
  if (abs(sam_table(t,"Government","Diff"))>zero_tol,
    display "Government account check      FAILED";
  else
    display "Government account check      PASSED";);
  if (abs(sam_table(t,"Rest of World","Diff"))>zero_tol,
    display "Rest of World account check   FAILED";
  else
    display "Rest of World account check   PASSED";);
  if (abs(k_check(t))>zero_tol,
    display "Capital account check         FAILED";
  else
    display "Capital account check         PASSED";);
  display   "------------------------------------";
);
display "----------------------------------------";



* if there is no policy file or baseline file then there is no cv to calculate
$if %policy_file%==0 $goto extra_code_start
$if %gdx_baseline_file%==0 $goto extra_code_start



* ------------------------------------------------------------------------------
* solve for equivalent variation
* ------------------------------------------------------------------------------

* solves for the willingness to pay for to forgo the policy change measured as
* the change in expenditures required to achieve the policy case welfare level
* at the baseline prices

parameters
  w_policy(r,h)         welfare in policy simulation
  exp_diff(t,r,h)       expenditure difference between baseline and ev case
  ev(r,h)               equivalent variation - finite horizon
  ev_annual(r,h)        equivalent variation - annual step and finite horizon
  ev_inf(r,h)           equivalent variation - annual step and infinite horizon
  c_ev(t,r,h)           equivalent variation consumption path
  leis_ev(t,r,h)        equivalent variation leisure path
  step_multiplier       multiplier for annual linear interpolation;

variables
  lambda_t(r,h)         terminal shadow price on consumption

equations
  aux_ev(r,h)           auxilary - equivalent variation constraint
  aux_lambda(t,r,h)     auxilary - lambda definition based on terminal value;

* save policy welfare to be used later in solving for ev
w_policy(r,h) = w.l(r,h);

* get the baseline prices
execute_load "%gdx_baseline_file%" pa, pk, pr, pl, pfx, pkt, kt, pr_ex, incadj,
                                   pres, res_ex, res, cpi, pr_ex_agg, tc;

* fix the prices at the baseline levels
pa.fx(t,r,s)      = pa.l(t,r,s);
pk.fx(t,r)        = pk.l(t,r);
pr.fx(t,r)        = pr.l(t,r);
pl.fx(t,r)        = pl.l(t,r);
pfx.fx(t)         = pfx.l(t);
pkt.fx(r)         = pkt.l(r);
kt.fx(r)          = kt.l(r);
pr_ex.fx(t,r,s)   = pr_ex.l(t,r,s);
incadj.fx(t)      = incadj.l(t);
pres.fx(t,r,s)    = pres.l(t,r,s);
res_ex.fx(t,r,s)  = res_ex.l(t,r,s);
res.fx(t,r,s)     = res.l(t,r,s);
cpi.fx(t)         = cpi.l(t);
pr_ex_agg.fx(t,r) = pr_ex_agg.l(t,r);

* constraint for ev calculation is that welfare must equal policy level
aux_ev(r,h)..
  w_policy(r,h) =e= w(r,h);

* calculate lambda based on the terminal value that equates welfare
aux_lambda(t,r,h)..
  lambda(t,r,h) =e= (1/(1+rho(r,h)))**(card(t)-ord(t))*lambda_t(r,h);

* define the model for computing ev
model calc_ev /
  zp_c.c
  zp_cm.cm
  zp_cme.cme
  zp_cene.cene
  zp_cen.cen
  zp_cl.pcl

  hde_cm.pcm
  hde_cme.pcme
  hde_cene.pcene
  hde_cen.pcen
  hde_c.pc
  hde_cd_trn.cd
  hde_cd_en.cd
  hde_cd_m.cd
  hde_cd_ele.cd
  hde_leis.leis

  foc_cl.cl

  uce_c.uc_c
  uce_cen.uc_cen
  uce_cme.uc_cme
  uce_cm.uc_cm
  uce_cl.uc_cl
  uce_cene.uc_cene

  aux_w.w
  aux_ev.lambda_t
  aux_lambda.lambda
/;

* iteration limit
calc_ev.iterlim = 100000;

* use model scaling
calc_ev.scaleopt = 1;

* use the path options file
calc_ev.optfile = 1;

* solve the mcp problem
solve calc_ev using mcp;

* save the equivalent variation consumption and leisure paths
c_ev(t,r,h) = c.l(t,r,h);
leis_ev(t,r,h) = leis.l(t,r,h);

* get the baseline consumption and leisure paths
execute_load "%gdx_baseline_file%" c, leis;

* calculate the present value reduction in expenditures required to reach
* the policy case welfare level at the baseline prices
exp_diff(t,r,h) = (c.l(t,r,h)-c_ev(t,r,h))*c0(r,h)*pc.l(t,r,h)
                  +(leis.l(t,r,h)-leis_ev(t,r,h))*leis0(r,h)*pl.l(t,r);

* calculate ev as the sum of present value reduction in expenditures
ev(r,h) = sum(t, exp_diff(t,r,h));

* cumulative sum -> sum_{i=0}^{interval}{i-1} used in the linear
* interpolation to annual time steps
set counter / 1*100 /;
step_multiplier = sum(counter$(counter.val le interval), counter.val-1);

* calculate ev for the model's finite time horizon but at an annual time step by
* linearly interpolating consumption and leisure paths between simulated years
ev_annual(r,h) = sum(t$(ord(t) ne card(t)),
                       interval*exp_diff(t,r,h)
                       +step_multiplier*(exp_diff(t+1,r,h)
                                         -exp_diff(t,r,h))/interval)
                 +sum(t$(ord(t) eq card(t)), exp_diff(t,r,h));

* calculate equivalent variation for an infinite time horizon assuming
* that the economy is in steady state
ev_inf(r,h) = ev_annual(r,h)
              +sum(t$tlast(t), exp_diff(t,r,h))
                *(1/(1-(1+(1+gamma)**(1/interval)-1
                         +(1+omega)**(1/interval)-1)
                       /(1+(1+rbar)**(1/interval)-1))-1)$(card(t) gt 1);

* output ev to results file
results.ap = 1;
put results;
loop((t,r,h), put "exp_diff",  t.tl, r.tl, "", h.tl, "", exp_diff(t,r,h) / );
loop((r,h),   put "ev",          "", r.tl, "", h.tl, "", ev(r,h)         / );
loop((r,h),   put "ev_annual",   "", r.tl, "", h.tl, "", ev_annual(r,h)  / );
loop((r,h),   put "ev_inf",      "", r.tl, "", h.tl, "", ev_inf(r,h)     / );
putclose;



* ------------------------------------------------------------------------------
* run extra code if provided
* ------------------------------------------------------------------------------

$label extra_code_start
$if %extra_code_file%==0 $exit

$include %extra_code_file%
